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FOREWORD 


This document was prepared by F. C. Kopper, Dr. G.J. Sturgess and Dr. P. Datta 
of Commercial Engineering, Commercial Engine Business of Pratt & Whitney, 
United Technologies Corporation, East Hartford, CT. It describes a study to 
develop and verify computational methods for predicting local heat transfer 
and coolant pressure drop within the coolant passages of rotating turbine 
blades. This work was carried out under the sponsorship of the National 
Aeronautics and Space Administration under Contract NAS3-23691 . This report 
describes a part of Phase I, Task III of this contract. The work was performed 
under the direction of Dr. F. Yeh, NASA Project Manager, and Mr. F. Kopper, 
Pratt & Whitney Program Manager. 
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1.0 SUMMARY 


The National Aeronautics and Space Administration (NASA) is sponsoring an ex- 
perimental and analytical program to investigate heat transfer within rotating 
turbine blades. The program is to generate flow, heat transfer, and pressure 
drop data for rotating passages simulating blade passages. The passages are 
multipass in configuration, with and without turbulators on the walls, to 
simulate the configurations and conditions expected in the first stage blades 
of advanced aircraft gas turbines. The experimental results are to be manipu- 
lated to provide data correlations for an empirically-based design system. The 
feasibility of developing computational fluid dynamic techniques for calculat- 
ing the flows in such passages is to be explored. 

This report presents a part of the planned Phase I, Task III effort which 
covers analysis and data correlation. 

The Pratt & Whitney 3D-TEACH (Teaching Elliptic Axisymmetric Characteristics 
Heuristically ) CFD (Computational Fluid Dynamics) code was selected as a suit- 
able vehicle for modification to meet the needs of the program. This state-of- 
the-art viscous flow computational fluid dynamics computer code has been 
revised to account for rotating internal flows. The modifications made have 
been evaluated for flows characteristic of those expected in the application. 
For this purpose, experimental studies extant in the literature were used. The 
results, while encouraging, reveal some limitations for the present codes. 



2.0 INTRODUCTION 


In current design technology there is a need to predict the effects of rotation 
on local heat transfer and pressure loss within coolant passages of turbine 
blades, particularly blades that have multipass coolant passages, such as shown 
in Figure 2.1. Accurate calculation of these quantities Is very important as 
it can lead to an improvement in the reliability of predicting turbine airfoil 
temperatures and ultimately, blade life. Accurate estimates of temperature and 
life permit available cooling air to be used most effectively and reduces its 
negative impact on turbine performance. 
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Figure 2.1 Blade Internal Geometry (JT8D) 


In recognition of these requirements, the National Aeronautics and Space 
Administration (NASA), through its Lewis Research Center and under the aegis 
of the HOST program, is sponsoring an experimental and analytical study of 
heat transfer and pressure loss in rotating multipass passages with configura- 
tions and dimensions typical of modern turbine blades. 

Computational fluid dynamics techniques, although far from completely devel- 
oped, have shown potential for calculating the Internal flows in the gas 
turbine engine. The complexity of the flow within turbine blades suggests that 
comprehensive empirical correlations of data will be difficult to obtain. 
Therefore, there is strong motivation to apply computational techniques to a) 
understand the flow development, b) assist in determining suitable correlations 
for the experimental data, and c) eventually, provide a quantitatively-accurate 
calculation procedure. 
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2.1 Objectives 


The overall objective of the NASA program is to develop and verify improved 
analysis methods that will form the basis for a design system which will result 
in efficient turbine components with improved durability. This objective is to 
be achieved through three program elements. The first element is to establish 
a comprehensive experimental database that can form the basis of an empirical 
design system. The second is to develop computational fluid dynamic techniques 
for this application. Finally, the third element is to analyze the information 
in the database with mathematical modeling to devise a suitable design and 
analysis procedure. The approach is illustrated in Figure 2.2. This report is 
concerned with part of the second of these elements - the computational fluid 
dynamic aspects. 





Figure 2.2 Program Structure 


The specific objectives of the work presently reported were to: 

1. Select a baseline CFD computer code 

2. Assess the limitations of the baseline code 

3. Modify the baseline code for rotational effects 

4. Verify the modified code against benchmark experiments in the literature 

5. Identify shortcomings in the code as revealed by the verification. 
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3.0 COMPUTATIONAL APPROACH 


Accurately predicting the local coolant-side heat transfer coefficients, 
coolant temperature rise and pressure drop over a wide range of operating 
conditions and geometries in the cooling passages Is a formidable task because 
of the turbulent, three-dimensional, elliptic nature of the flow. Further 
complications are introduced by the effects of rotation and turbulence - 
promoting devices frequently introduced into the passages. Figure 2.1 
illustrates these features in a typical passage. 

3.1 Basic Procedure 

Computational fluid dynamics methods in general fall Into four main classes, 
namely : 

1. Incompressible potential flow solutions 

2. Solutions based on viscous-inviscid interaction approaches 

3. Parabolic time-marching finite difference methods for viscous flows 

4. Elliptic Iterative finite difference methods for viscous flows 

It can be appreciated that this list of classes is one of increasing physical 
realism and flexibility to handle more complex problems. It also is a list 
representing increasing difficulty, complexity, and cost. The character of the 
flow field described above, and to be calculated, dictates an approach from 
Class 4. 

Within Class 4, the TEACH-code was selected for the present program because it 
represented the only Class 4, generalized, three-dimensional, viscous, elliptic 
flow code currently available at Pratt & Whitney that was sufficiently devel- 
oped to be considered. The approach used in this code for turbulence management 
(time -mean /wholly statistical) represents a technically reasonable and com- 
putationally economic approach in an engineering environment (Ref. 1). 

Th'e computational approach selected is one being developed by Pratt & Whitney 
for three-dimensional, viscous, reacting, elliptic flow calculations for the 
combustor and other gas turbine applications. The procedure is known by the 
acronym 3D-TEACH. It is one of a family of such codes being developed by Pratt 
& Whitney, and its generic solution approach Is well-known and accepted in the 
industry. The acronym TEACH (Teaching Elliptic Axi symmetric Characteristics 
Heuristically ) represents a generic solution technique (Ref. 2) and these codes 
represent current production state-of-the-art calculations in terms of 
equations solved, physical models used, discretization of the equations, and 
the solution algorithms. They are not perfect, but are a marked advance from 
one-dimensional flow calculations and global modeling that formed the previous 
capability. The structure of the codes has been made such that modular 
replacement can be carried out as better models and solution algorithms are 
developed, while the basic framework and operational features remain. 
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The Pratt & Whitney 3D-TEACH code is a generalized aerothermal , fluid dynamic 
solver for steady, three-dimensional, elliptic, turbulent, reacting flows. The 
approach stays within the framework of continuum mechanics and uses a statis- 
tical description of turbulence, coupled with the accepted Eulerian description 
provided by the Navier-Stokes equations of motion. Closure to the resulting 
time-mean equations is provided by turbulence modeling of the eddy viscosity 
type. The modeled partial differential equations are manipulated into a general 
form that permits a single solution algorithm to be used for a numerical pro- 
cedure. A hybrid (upwind/central) finite differencing scheme is used to dis- 
cretize the equations. An outline of the solution procedure is given in 
Appendix Al. In original form, the code was derived at Imperial College, 

London, England, as a teaching aid. 

3.2 Baseline Code and Its Features 

The Pratt & Whitney 3D-TEACH code is being developed for application to the 
combustion chamber of the gas turbine engine (Ref. 3). In cylindrical coordi- 
nates, it solves the general steady flow equation. 
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where overbars denote time-averaged quantities, and. 


(3.1) 


p = gas density 


u, V, w 


gas velocities in x, r, © directions, respectively 


i = any of the independent variables 

r eff,e = an appropriate turbulent exchange coefficient, depending on what 

t represents 


= a so-called "source term" which lumps together all other items 
in a given equation not included in the first six terms of 
Equation 3.1. 


This general form was adopted so that a single solution algorithm could be 
used in the numerical procedure, which is given in Appendix Al. By way of 
example. Figure 3.1 gives the values of some of the parameters in Equation 3.1. 
Those shown are relevant for the present application. However, the code also 
contains information relevant to its combustion application, including, for 
example, species equations with source terms to describe the reaction rate of 
a fuel, etc. This extraneous information was stripped from the code to estab- 
lish a baseline code for use in this turbine study. 
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Figure 3.1 Summary of Some of the Equations Solved in 3D-TEACH 
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In Figure 3.1, the terms K and e respectively represent the specific kinetic 
energy of turbulence and its rate of dissipation. Thus, the code solves trans- 
port equations for K and e. These arise from the turbulence model incorporated, 
which is the so-called two-equation, or K-e, model. The turbulence model pro- 
vides expressions for the Reynolds stresses in terms of calculable quantities. 

Reynolds stresses can consist of two parts - a shear stress and a normal 
stress. The normal stress is obtained simply from the fluctuating dynamic 
pressures, while Boussinesq's analogy is used to relate the shear stress to 
the velocity gradient through an eddy viscosity ut- Thus for incompressible 
but variable density flow, the Reynolds stress can be expressed in compact 
tensor notation, as: 
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where primes represent randomly fluctuating values, and. 


<5-jj = Kronecker delta 

K = 1/2 ( q '2 + y '2 + - 2 ) 


(3.2) 


(3.3) 


The eddy viscosity Mt is obtained by dimensional arguments from the Prandtl- 
Kolmogorov definition. 


where Cm is a constant of proportionality, equal to 0.09. 

As the flow field is found from an effective turbulent eddy viscosity, it is 
convenient to also base the turbulent heat transfer on an effective thermal 
diffusivity. 

Eddy diffusivity gives for the flux of a scalar. 


where 

e = a scalar quantity, i.e. temperature 
r t = turbulent eddy diffusivity. 


The eddy diffusivity is found from the ratio of turbulent kinematic eddy vis- 
cosity v£, to a turbulent Prandtl number j5 t , that is specified as input.' 
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For the eddy viscosity approach to turbulence modeling, the two-equation ap- 
proach is the most general and sophisticated representation, and it is not 
computationally expensive. The sophistication comes from the use of differen- 
tial equations to describe both the velocity scale and length scale to which 
eddy viscosity is assumed proportional, rather than relying on an a priori 
scale specification as in the mixing-length approach. 

There are a number of limitations to the model. The K-equation is exact, but 
the modeling used in producing the dissipation equation is known to be shaky. 
More Importantly, use of the gradient diffusion idea itself has been long 
challenged. There are objections to the assumption that the Reynolds stresses 
depend on just the local mean rates of strain, as well as to the assumption 
that the stresses are proportional to these rates of strain. The "constant of 
proportionality" C really depends on the ratio of local production and dissi- 
pation of turbulence energy, and this ratio is not actually a constant. A 
further weakness is the adoption of a single velocity scale at a point in the 
flow, although this scale can vary from point to point. The implication of a 
single scale is that the turbulence Is isotropic. Real turbulence usually has 
some degree of anisotropy. Certain flows result in turbulence which is highly 
anisotropic; for example, flows that are swirling or which have strong curva- 
ture in the streamwise direction. The velocity and length scales have to be 
the same order of magnitude as the mean field motion. This is only true for 
flows dominated by simple shear forces; buoyancy forces, for example, have 
their own separate scales. It is implied that the turbulent motions have a 
small scale compared to that over which the magnitude of a diffusing quantity 
changes significantly. Most of the larger eddies in a turbulent flow do not 
satisfy this condition, whether the eddies are coherent or not. Thus, material 
can be transported by vortical motion against the gradient of temperature. 
Another "action-at-di stance" which is removed from consideration by relying on 
local mean rates of strain is the effect of "flow history" on turbulence 
structure. 

The shortcomings of the turbulence model and their relevance to the internal 
flows of the turbine blade will be returned to later. Two other features of 
the baseline code have to be described because they are of direct relevance to 
the surface heat transfer calculations ultimately to be made in the present 
application. 

The first feature is related to the treatment of curvilinear surfaces, and it 
arises from the use of finite differences to approximate the partial differen- 
tial equations to be solved. 

The finite difference analogue of the differential equations is obtained by 
overlaying a computational mesh on the flow domain to be calculated, and 
obtaining the basic finite difference form of the partial derivatives for every 
node of the mesh from a control volume approach. The finite difference expres- 
sions, when substituted back into the differential equations, yield a set of 
linearized, algebraic equations (Appendix A1 ) for every node of the mesh. In 
finite difference methods, the mesh is constrained by the coordinate system 
used; therefore, either the coordinate system is chosen to fit the geometry, 
or the geometry is "discretized" to fit the coordinate system. 
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The 3D-TEACH code is a generalized code. Therefore, the coordinate system is 
fixed, and the boundary geometry is made to conform. All curvilinear surfaces, 
therefore, have to be represented through "stair-steps". 

The use of stair-step geometries has a number of implications. First, surface 
areas are not correct. Thus, Irrespective of physical modeling and numerical 
accuracy, calculation of wall shear stress and surface heat transfer on sur- 
faces not aligned with the mesh can never be correct. Second, adequate repre- 
sentations of the geometry bounding the flow to be calculated usually require 
more computer storage than is available on the current generation of computers. 
Mesh refining to control numerical diffusion (to be discussed) is not there- 
fore possible, and furthermore, in some circumstances the calculated flowfleld 
may be influenced by the geometric representation. 

In order to calculate surface heat transfer, some treatment of the wall bound- 
ary layers is necessary. The K-| turbulence model, described above, involves 
many assumptions and has known shortcomings; In addition, it is essentially 
for high Reynolds numbers. It is appropriate to examine the suitability of 
this model for the boundary layer calculations needed. 

With the no-slip boundary condition established at Impervious walls, the local 
velocity adjacent to solid surfaces becomes low. As walls are approached there- 
fore, the local Reynolds number of the flow based on local velocity and dis- 
tance from the wall Is very small, and viscous stresses become significant. 

The large-scale turbulence structure is likely to be directly influenced by 
viscosity, and the presence of the solid boundary will distort the turbulence 
structure. 

Near to a solid surface, the turbulent eddies are of a size comparable with 
the distance from the surface since they cannot be larger than this distance. 
This physical situation is that in the regions of the flow where the turbulence 
is produced. In the high Reynolds number production regions, the turbulence- 
producing eddies are widely separated in wave-number space from regions where 
this energy is dissipated. In the wall region, the dissipative motions are 
Influenced directly by the mean strain rate, so that the fine scale eddies 
share the anisotropy of the large scale eddies and are no longer isotropic. 
Characterization of such turbulence in terms of a single length and time scale 
is a questionable assumption. Because the Reynolds stresses fall to zero on 
the solid surface, very large gradients of turbulence properties exist normal 
to the wall. These gradients cause significant turbulent transport of energy 
toward the wall, so that the inner layer of the boundary layer is not in local 
equll ibrium. 

The existence of simple boundary conditions leads to simply geometries and 
relatively simple mathematics for boundary layers. However, the turbulence of 
boundary layers is just as complicated, or more so, as that of high Reynolds 
number flow. The two-equation turbulence model is inappropriate for boundary 
layers because it neglects any direct viscous effects on the turbulence struc- 
ture. Most of the assumptions inherent in the two-equation model are violated, 
i.e., single length and time scales, local equilibrium, and homogeneity. 
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To take account of the direct influence of viscosity a calculation has to con- 
tinue solving the conservation equations right through the viscous sub-layer 
of the boundary layer, where all the viscous effects on turbulence dominate. 

If the two-equation model of turbulence will not do, the physics could be 
accounted for by an appropriate higher order closure. However, higher order 
closures are inherently computationally expensive to use, and become almost 
prohibitively expensive in the case of boundary layer calculations. Because of 
the severely steep gradients normal to the surface in boundary layers, not 
less than 15 cross-stream grid nodes are required for the sub-layer and a 
total of 100 may be required for the complete layer. 

The TEACH code was originally conceived for the calculation of complete 
fields, and not just boundary layer flows. For the reasons given above, an 
attempt to simultaneously calculate the bulk flow field and the boundary layer 
flows would result in an enormous computer storage requirement, and would be 
extremely expensive to perform. To circumvent these difficulties, the TEACH 
codes use special treatment at the walls. 

It was argued that the main flow field is not Influenced to first order by the 
details of the flow at the walls. The profile "universality" feature of bound- 
ary layers was then invoked to provide an easy description of the flow near 
solid boundaries. A "law of the wall" approach was used to conveniently link 
the wall to the near-wall nodes of the finite difference solution grid. This 
is done very inexpensively. The approach is described in Appendix A2. 

The wall -function method is a standard procedure in all the TEACH codes, and 
is satisfactory for most purposes. However, deviations from "universal" laws 
do arise, and the wall-function approach is not valid for separated flow 
regions, flows with large density gradients, boundary layers with strong cur- 
vature, and most Importantly, large pressure gradients in the direction of 
flow. 

It should be appreciated that the TEACH code Is not analytically exact but is 
an approximate solution procedure. Approximations are introduced at several 
stages. However, of particular importance is the use of finite differencing to 
discretize the equations to be solved. 

The accuracy of a differencing scheme can be judged from the order of the 
terms of an equivalent Taylor series that have been retained in the expansion. 
Unfortunately, the requirements of numerical stability are opposite to those 
of accuracy with respect to these terms. The spatial differencing of the con- 
vective terms of the conservation equations in an Eularian coordinate system 
can result in numerical diffusion occurring. Use of a higher order differencing 
scheme eliminates or significantly reduces this diffusion. However, use of 
central differencing, for example, introduces an oscillatory behavior into the 
solution. This "wiggling" can lead to nonphysical behavior (Ref. 4). The use 
of an upwind or donor-cell technique eliminates wiggling; however, this is 
accomplished by the Introduction of a diffusive-like term into the difference 
equations. Thus, while "numerical damping" suppresses oscillation, it leads to 
significant additional diffusion of the convected parameter. Therefore, a 
severe restriction can be placed on the quality of quantitative prediction 
(Ref. 5). 
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It can be argued for use of upwind differencing in regions where convection 
strongly dominates streamwise diffusion, that the local upstream values of the 
field variables are swept downstream virtually unchanged, whereas in high- 
diffusion regions the form of the relatively small convection terms is not 
important. In regions where the two transport mechanisms are comparable, a 
switch to more accurate central differencing for convection or use of a suit- 
ably weighted combination of central and upstream differencing can be used. 

This somewhat narrow view of complex flows has led to the appearance and use 
of a popular and successful hybrid central /upwind differencing scheme (Ref. 6). 
This scheme is currently used in the TEACH codes described in Appendix A1 . The 
scheme uses central differencing for convection and diffusion fluxes when the 
absolute value of the Peclet number for the control volumes existing about 
grid nodes is less than, or equal to, two; upwind differencing for convection 
fluxes and neglect of diffusion fluxes is used otherwise. Peclet number defines 
the relative importance of convective and diffusive transport. 

To successfully use the hybrid differencing scheme for complicated flows, care 
must be taken in establishing the computational grid upon which the calcula- 
tions are performed. The approximations of the algebraic expressions used to 
represent the partial differential equations becomes asymptotically exact as 
the distance between the nodes set up by the grid, and used to link the 
algebraic expressions, is reduced. In the limit, the number of nodes can be 
increased until an asymptote to the solution to the differential equations is 
achieved. In practice, this increase is limited by computer storage and the 
cost of the calculation. However, it is not just the number of nodes that are 
used which is important in determining the accuracy of a solution, but also 
the distribution of those nodes within the flowfield to be determined (Refs. 7 
and 8). This nodal distribution is Important because whenever curvature of the 
flow in the streamwise coordinate direction exists, a truncation error arises 
In the solution (Ref. 9). There is also a problem in multidimensional flows of 
streaml ine-to-gri d skewness (Ref. 10). With upwind differencing, these effects 
start to have a damaging effect on solution accuracy when the Peclet number 
exceeds two. 

It has been concluded (Ref. 3) that the hybrid finite differencing scheme, 
although yielding physically realistic solutions in all circumstances, intro- 
duces excessive numerical diffusion for many two-dimensional flows, and for 
all three-dimensional flows. This is because present computer storage is not 
sufficient to permit local adjustment of the grids as described above, except 
for the simplest of flows. Thus, the solution accuracy is presently controlled 
by the numerics, rather than the hierarchy of physical modeling. 

3.3 Suitability of Baseline Code 

The salient features of the baseline code have been described. Some limitations 
and shortcomings were highlighted. The suitability of the baseline code for 
modification to a form for the intended application is now discussed. (The 
modifications to be made are illustrated in Figure 3.2.) 
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Figure 3.2 Modification Procedure to Develop Code for Turbine Blade Analysis 


There are four features of the baseline code that give rise to concern for 
application to turbine blade Internal flows: 

1. Suitability of the turbulence model for calculation of the bulk flow 
f i el d ; 

2. Stair-step treatment of curvilinear surfaces; 

3. Use of wall functions to represent boundary layers; 

4. Effects of numerical diffusion. 

3.3.1 Turbulence Model 

It has been found experimentally that general streamline curvature in the 
plane of the mean shear produces large changes In the turbulence structure of 
the shear layers. These effects can occur not only in boundary layers but in 
any shear layer where the streamlines have a component of curvature in the 
plane of the mean shear. Turbulence in the boundary layer on a highly convex 
surface may be nearly eliminated, while on highly-concave surfaces momentum 
transfer by quasi -steady longitudinal vortices dominates the ordinary turbu- 
lence processes. Thus, the changes in turbulence are not only quantitative but 
also qualitative. The changes are usually an order of magnitude more important 
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than normal pressure gradients and other explicit terms appearing in the time- 
mean equations of motion; this is not so for flows dominated by secondary 
flows. This effects on local Reynolds stresses do not appear as soon as the 
curvature is imposed. This means flow "history" effects are important. 

The streamline curvature is a distortion that gives rise to extra rates of 
strain in addition to that due to simple shear. The Reynolds stresses are 
changed by a factor (1-as) where 5 < a < 15, and 

s = Q/r/( a-/ar) 

Now, consider the flows sketched below where in the first case the angular 
momentum increases outward from the flow center of curvature, and in the second 
case where it decreases. 

I II 




Angular momentum Angular momentum 

increasing with radius decreasing with radius 

It can be appreciated from the equations above that for the first case the 
Reynolds stress is reduced, the turbulence is suppressed, and the flow curva- 
ture is stabilizing. For the second case, the turbulent motion is strongly 
augmented, and the large eddies develop into longitudinal vortices. 

The two-equation turbulence model with its assumption that the one-point, two- 
variable correlations for Reynolds stresses and turbulent scalar fluxes are 
directly proportional to th mean gradients implies that the eddy viscosity is 
isotropic. The model therefore displays equal sensitivity to both primary and 
secondary strains; certainly it also has no way to account for flow "history." 
The flow behavior is not well-predicted by such mean-field closure methods. 

The flow through a blade cooling passage (see Figure 2.1) will involve at least 
one 180-degree sharp turn, and more with multi-pass configurations. Flows in 
ducts with curvature produce secondary flow motions Induced by the transverse 
pressure gradient due to lateral curvature of the main flow. In addition, 
secondary motions arise for turbulent flows in non-circular ducts due to 
Reynolds stress gradients in the cross-stream direction, particularly for thick 
inlet boundary layers. As the flow proceeds through the passage, the effects 
become more pronounced. The flow curvature produced by the turns will also 
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result in the stabilizing and destabilizing changes in turbulence structure 
described above. However, in this case the turns are likely to be so sharp 
that the cross-stream flow is probably dominated by mean pressure forces 

(Ref. 11), so that the effects of the turbulence model on mean profiles should 
be smal 1 . 


For the model to bemused in the present experiments (Figure 3.3), the passage 
high aspect ratios in the flow direction and the modest curvatures used suggest 
that boundary layer effects could be quite Important. If this is so, the code 
could not be expected to do well with the isotropic K-e turbulence model. 


COVER 



Figure 3.3 Coolant Passage Heat Transfer Model 
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3.3.2 Treatment of Curvilinear Surfaces 


The stair-step representation of the modest curvatures shown for the model in 
Figure 3.3 will introduce a much larger momentum deficit into the flow than 
would exist physically. If the flow in the model were to separate in leaving 
the turns, this deficit would not be so critical. However, the turns are so 
modest that significant flow separation on the inner surface of the turn is 
unlikely to take place. This being so, the boundary layers in the legs away 
from the turns will be calculated incorrectly from those actually existing. 

The flow in an actual blade (see Figure 2.1) is likely to be better calculated 
than that in the present model. This is because the baffle in a real blade is 
much thinner than in the model, such that a rectilinear representation of it 
does not introduce any significant momentum deficit into the flow. (The turning 
vanes suppress separation.) 

3.3.3 Use of Mall Functions to Represent Boundary Layers 

The flow through the model has large density gradients, significant pressure 
gradients in the direction of flow, boundary layers with moderate to strong 
curvature, and possibly small regions of separated flow. These are all situa- 
tions in which the wall function approach could be expected to fail. However, 
the significance of such failures is likely to be overshadowed by the limita- 
tions of anisotropic turbulence, the crude representation of the bend curva- 
tures, and the effects of numerical diffusion. 

3.3.4 Numeri cal Pi ff usi on 

The calculations to be made are elliptic in character. Therefore, the imposed 
downstream boundary conditions have to be sufficiently removed from the region 
of interest such that they do not influence the calculation in this region. 

For the model shown in Figure 3.3, this means that to obtain good calculations 
of the flow in the first two legs of the passage, three legs must be calculat- 
ed, and so on. Thus, the calculation domain becomes extensive. An extensive 
calculation domain can demand a large number of grid lines to cover it with an 
acceptable degree of resolution. 

The maximum array size currently available on Pratt & Whitney's computers is 
of the order of 30x40x40, or 48,000 nodes. A considerable number of grid lines 
are consumed in representing the turns with stair-steps and in defining the 
flow passages. The three-dimensional grid necessary to adequately represent 
and fully cover a complex and extensive calculation domain quickly consumes 
the available computer storage. Thus, there is little opportunity to refine 
the mesh to minimize numerical diffusion, as is possible in two-dimensional 
calculations (Ref. 12). Therefore, considerable numerical diffusion will in- 
evitably be present in the calculations. This is a problem with all three- 
dimensional viscous flow calculations made on the current generation of general 
purpose computers (Refs. 3 and 13). 

3.3.5 Summary 

The selected baseline code has the basic features essential for the stated 
objectives. It is therefore a suitable vehicle for modification for rotational 
effects. However, thehe are limitations and shortcoming to the code that will 


impact the quantitative accuracy achievable. Fortunately, the structure of the 
code is modular so that some of the deficiencies may easily be remedied as 
better models come along. This applies specifically to the turbulence model 
and the finite differencing scheme. Also, it is anticipated that computer 
storage will soon cease to be a problem (although the cost of solution will 
then need to be addressed by speeding up the calculation procedure). The pro- 
blem of a three-dimensional curvilinear representation is much more difficult 
and a near-term answer is not at hand. 5 

3.4 Modifications for Rotation 

The equations to be solved are the three momentum equations, the continuity 
equation, the energy equation, an equation of state, and the transport equa- 
tions for K and e. The modifications to be made involve Incorporating rota- 
tional effects. If the blade is placed at radius R and is rotated with steady 
angular velocity additional inertial accelerations will arise. These are 
due to Coriolis and centrifugal forces, 2^27 and -g&R respecti vely, 
where V is the bulk flow velocity in the cooling passage, and where the 
overbar denotes a time-mean quantity as before and denotes a vector quantity. 

The coolant passages, shown in Figure 3.3, are most conveniently represented 
by a Cartesian coordinate system x, y, z, fixed on the "blade." The additional 
terms due to rotation are most conveniently added to the general equation 
(Appendix Al, rewritten in Cartesian form) as additional source terms, n. 
Therefore, the Coriolis and centrifugal forces must be derived in terms of the 
rotating coordinate system shown in Figure 3.4. 



Figure 3.4 Coordinate System for a Rotating Frame of Reference 
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Coriolis forces can contribute directly to the development of secondary flows 
and indirectly to affect the turbulence structure. The rotation has a stabili- 
zing-destabilizing effect on turbulence, similar to that described for bulk 
flow curvatures described above. It has been observed experimentally (Ref. 14) 
in rotating channel flow that on the trailing (suction) side rotation suppres- 
sed turbulence production and on the leading (pressure) side rotation caused 
the appearance of Taylor-Gorl ter-type vortices. 

There is no contribution to the continuity equation. The contribution to the 
momentum equations is the sum of the Coriolis and centrifugal terms. The 
velocities J2 and 7 are resolved into their components, (ifi x + jfty + kS2 z ) 
and (iu + jv + kw) respectively, to give the additional sources 50, n . Unit 
vectors in the major coordinate directions are represented by i, j and k. The 
momentum equation is shown in Appendix A3, while the sources are given in 
Table I. 

TABLE I.- ADDITIONAL SOURCE TERMS DUE TO ROTATION 
Equation t 

x-momentum u p ( v-z J2 X + x z ) z - p (w - xfl y + y fi x ) n y 

y-momentum v p(w-xS2 y + yJ2 x )S2 x -p(u - yfi z + z J2 y ) fi z 

-pu-n z -pv?n x 

w p (u - y fi 2 + z n y ) S2 - p (V - z fi x + x n z ) J 2 X 

+ p uJ2 y - p Vfi x 

h p u [J2 y (n^ -n y x) +n 2 (n x z - n z x)]| 

-|pv[s2 x (n y x-n x y) + n z ( n y z - ^1]} 

- 1 ? w [ n x ( n 2 x - n x z) + fi y ( 


Similarly, there is a stagnation enthalpy term p^ ^ R that gives rise 
to a source term in the energy equation, as Table I shows. . 

No additional terms arise in the K-equation. An assumption of the turbulence 
model, that the turbulence Is in equilibrium so that production and dissipation 
of turbulence kinetic energy are in balance, is invoked to infer that there 
are no additional terms in the e-equation either. This means that although 
secondary flow development due to Coriolis forces should be calculated through 
the modified equations of motion, the effects on the turbulence structure will 
not be calculated. 


z-momentum 


enthal py 
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The passage angle e In Figure 3.4 Is accounted for in the resolving of ft into 
Its components. The position vector r in Figure 3.4 is automatically taken 
care of by working in terms of the angular velocity R. 

Details of the derivations are given in Appendix A3. Buoyancy effects are not 
accounted for as they should be small in this application. 

3.5 Exploratory Calculations 

Following debugging of the programming changes to Incorporate the rotational 
terms shown in Table I, it was considered desirable to establish if the 
modifications gave a calculated flow behavior that was at least physically 
realistic. To this end, a simplified passage geometry was considered, as shown 
in Figure 3.5. 



Figure 3.5 Schematic of Coolant Passage for Exploratory Calculation 


In the simplified geometry a single outflow and inflow leg was considered. The 
grid used (33x11x5) is shown in Figure 3.6, where it is seen to be extremely 
coarse. However, this grid utilized the biggest virtual machine then available 
on Pratt and Whitney's computer system (4 Megabytes; currently 12 Megabytes). 
It illustrates the current difficulties associated with three-dimensional 
calculations. The calculations that result from this grid have to contain 
large amounts of numerical diffusion and are certainly not grid-independent. 
Although the flow resolution attainable with this grid is also poor, the 
calculations are considered adequate to ascertain whether or not the modified 
code is behaving as it should. 
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Figure 3.6 Grid for Exploratory Calculations 


The operating conditions for the exploratory calculations are shown in Table 


TABLE II.- OPERATING CONDITIONS FOR EXPLORATORY CALCULATIONS 


Passage cross-section 
Passage length 
Inlet air temperature 
Air pressure 
Coolant mass flow rate 
Wall temperature 
RPM 

Reynolds number 
Rossby number 


1 .27 x 1 .27 cm 
1 5.24 cm 
1 5.9°C 
10 bar 
0.0077 kg/s 
94°C 
600 
30,000 
0.174 


(0.5 x 0.5 inch) 
(6 inches) 

(60°F) 

(10 atmos) 

(0.017 lb/sec) 
(200°F) 

(basel ine) 


The table represents a "smooth-wall " test condition for the model of Figure 3.3. 
The passage is rotated about the y-axis. 


The results of the calculations are presented as flow visualizations by means of 
streaklines. Streaklines are the computational analogue of fine aluminum tracer 
particles momentarily illuminated by a finite-thickness sheet of laser-light as 
used, for example, in water tunnels. Streaklines are not streamlines, nor are 
they vectors, although they have some characteristic of both. The TEACH code 
can generate either random streaklines or uniform lines of origin; in the 
present case the random mode is used. 


Figure 3.7 shows the flow In the passage at cross-sections of 19.8 cm 
(7.8 inches) and 26.8 cm (10.56 inches) from the axis of rotation for a rota- 
tional speed of 600 RPM. In Figure 3.8, flows In the same sections are compared 
for 0, 60, 600 and 1900 RPM, and the flow through the complete passage length 
(x-y plane) is also shown. 


With reference to Figure 3.7, the development of secondary flows due to the 
Influence of Coriolis forces can be seen. In the outflow (away from the axis of 
rotation) leg of the passage, a pair of counter-rotating vortices develop. The 
vortex centers are shifted slightly from the centerline of the leg towards the 
pressure side, but the size and strengths of the vortices are about equal. As 
the vortex pair enters the turn, their rotational velocity is overcome and flow 
changes direction as mass is forced to the pressure side of the turn. The flow 
entering the passage inflow leg from the turn is therefore forced into a 
right-angled corner. As it escapes from the corner to begin flowing down the 
Inflow leg (toward the axis of rotation), the air has no choice other than to 
establish a single vortex. The action of the turn is thus to coalesce the vortex 
pair formed by Coriolis forces into a single vortex, completely filling the 
passage Inflow leg. Figure 3.8 demonstrates that increasing rotation 
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Figure 3.7 Test of Momentum Equation - Flow Visualization 


increases the strength of the single vortex in the inflow leg leaving the turn. 
However, the double vortex system begins to re-establish itself in the inflow 
leg as the flow proceeds toward the axis of rotation. The pair of vortices so- 
formed is Initially no longer symmetrical as the vortex originating from the 
coalesced pair from the outflow leg dominates. With rotation there is no evi- 
dence of flow separation in the inflow leg immediately following the turn. This 
is due to the vortex motion. 

These results are intuitively correct (Ref. 15) and show that the code is be- 
having soundly in a qualitative sense. Consistency checks were also carried out 
on this geometry to establish that the same results were obtained for all 
orientations. 
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4.0 EXPERIMENTAL DATA 


Although the modified code appeared to be performing in expected fashion as a 
result of the exploratory calculations, some form of quantitative verification 
was desirable before committing to extensive calculations for the test geometry 
shown in Figure 3.3. In principle, the verification testing should be conduct- 
ed for geometries as simple as possible, and against experiments that explore 
the overall relevant physical processes one at a time and not simultaneously. 
These two conditions are essential so that the effects of geometric modeling 
can be separated from the physical modeling, and so that shortcomings can be 
identified with a specific physical model. Ideally then, experiments of simple 
isothermal flow with rotation, stationary isothermal flow in sharp 180° bends, 
heat transfer in simple flow with rotation, and, heat transfer in stationary, 
sharp, 180° bends, would represent a desirable set of experiments against 
which to verify the performance of the modified three- dimensional code. 

4.1 Benchmark-Quality Experiments 

A series of experiments was defined above that would provide a suitable veri- 
fication. However, to be useful for this purpose, an experiment must also 
satisfy additional criteria. A qualified experiment that does satisfy these 
additional criteria is termed a benchmark -qual ity experiment. Benchmark -qual ity 
experiments are difficult to find. 

A benchmark -qual ity experiment is defined as follows: 

1. Minimum necessary flow dimensionality . Experiments in which the flows 
can be represented as one or two-dimensional are desired. Three- 
dimensional flow situations will be used only when specifically 
testing three-dimensional flow modeling capability. 

2. Well-behaved flows . Flows in which instabilities, periodicity, or 
changes in gross behavior occur as flow conditions change are avoided. 
For example, flows in which the location of reattachment points of 
separated flow regions could undergo significant shifts as Reynolds 
number is changed over the range of interest should be avoided 
(unless, of course, this is the flow feature being tested). 

3. Continuous variation of test parameters. Experiments should be 
conducted over a wide range of values of test parameters rather than 
at isolated sets of conditions. 

4. Known boundary conditions . Entrance and exit flow profiles must be 
specified as completely as possible. Velocity, temperature and 
pressure profiles are required. Concentration profiles are important 
for reacting flow experiments; initial droplet size, velocity and 
spatial distributions are required for two-phase flows. For assessing 
turbulence models, initial profiles of turbulence intensity and 
integral length scale are vital. 
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5. 


Progression In flow complexity . The ideal experiment for assessing 
flow models would consist of a series of experiments of increasing 
flow complexity such that the credibility of the analysis can be 
checked in stages. For example, the initial tests for a given flow 
geometry should be single-component, isothermal flow visualization 
tests which could be used to check the fluid mechanic aspects of the 
analysis. Two-component gaseous flows could be used to check the 
ability to predict mass diffusion; thermal diffusion could be checked 
using gases introduced with different initial conditions, reacting 
gaseous mixtures the next stages. A similar progression of experiments 
with two-phase flows can be constructed. 

6. Extensive 1 nstrumentation . Flow mapping experiments in which 

nonintrusive techniques are used to characterize flows throughout the 
chamber volume are highly desirable. Estimates of instrument precision 
should be made and possible sources of bias identified. Redundant 
measurements performed with different instruments are valuable. 

For the present purposes, extensive and complete verification is not the 
objective so that a full-range of experiments cannot be considered. Despite 
this, those experiments chosen should conform to the definition as closely as 
possible. Of course, it is most unlikely that any single experiment will con- 
form exactly to the ideal, so that some compromises have to be accepted. 

4.2 Selected Experiments 

An exhaustive literature survey was not made, and a total of only three 
experiments was selected. Two of these experiments explored rotational effects 
In isothermal flow and also covered the effects of passage aspect ratio 
(Refs. 17 and 18). The third experiment was concerned with heat transfer in a 
sharp, 180 bend, and had no rotation (Ref. 19). Although a more extensive 
comparison against experiments in the literature, and a literature survey 
would have been useful, such a comprehensive verification and survey were not 
a part of the present contract. 

4.2.1 Effects of Rotation 

Calculation of the effects of rotation on isothermal flow along a duct was 
investigated through the experiment of Moon (Ref. 17). Although carried out 
some years ago, it seems to have been a well-conducted experiment, and closely 
fits the description of a benchmark -quality experiment. 

The test-duct was rectangular, with a cross-section of 15.24 x 7.62 cm (6 x 3 
inches) and a length of 1.83 m (72 inches). Air was supplied uniformly at one 
end with a mean velocity of 16.5 m/s (54 ft/sec). The duct was arranged hori- 
zontally with the major axis of the cross-section vertical, and was rotated 
horizontally at a steady speed of 165 RPM about an axis passing through the 
center of the duct and perpendicular to the main flow direction. The arrange- 
ment is shown in Figure 4.1. The entering air was delivered through a contrac- 
tion and the supply duct contained honeycomb and screens to deliver a flat 
inlet velocity profile unaffected by rotation. The rotating test duct was 
contained within a cylindrical box to minimize disturbances. 
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Figure 4.1 Rotational Passage of Moon's Experiment 


Measurements were confined to the horizontal centerline at stations shown in 
Figure 4.1. The first station ( 1 M ) was located 0.46 m (1.5 feet) from the duct 
inlet, and the other stations were 0.30 m (1 foot) apart. A flattened hypoder- 
mic tube was used to take profiles of total pressure across the minor axis of 
the duct. Six static pressure taps were provided on each side-wall. Wall shear 
stress measurements were made by means of Preston tubes. Fluctuating veloci- 
ties and Reynolds stress profiles were measured with hot-wires. At the inlet, 
horizontal profiles of mean and fluctuating axial velocities were provided. 

4.2.2 Effects of Duct Aspect Ratio With Rotation 

Moon's experiment was conducted with a duct aspect ratio of 2:1. Moore 
(Ref. 18) used essentially the same apparatus to explore the effects of duct 
aspect ratio with rotation. This experiment was considered significant because 
passage aspect ratio is an important design parameter In practice. 

The modified duct section of the apparatus is shown in Figure 4.2. The spacing 
between the vertical side-plates now forming the duct was set to 1.9 cm (0.75 
inch) for passage aspect ratios (height/width) of 1:1, 4:1, and 7.33:1; while 
a value of 3.81 cm (1.5 inches) was used for an aspect ratio of 0.5:1. The 1.9 
cm (0.75 inch) duct width allowed fully-developed flow to be established. 
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Pratt and Whitney has sponsored a program of work at Arizona State University 
to provide data on pressure losses and heat transfer in sharply turning 
passages. Some of this material has been utilized for the present purposes 
(Ref. 19). 

Figure 4.3 shows the arrangement of the flow passage, and it can be seen that 
the 180° turn is really sharp. This is as it is in a real blade (see Figure 
2.1), and not as it is in the model to be tested as part of the present 
contract (see Figure 3.3). The turn also can involve either contractions 
(Wl/W2>1.0) or expansions (Wl/W2<1.0) of the passage. Passage aspect ratio can 
be varied as well, as can the end-wall clearance beyond the partition forming 
the sharp turn. The apparatus is stationary. The overall passage consists of 
an unheated plenum chamber provided with screens to ensure flow uniformity, 
the heated test-section, and an unheated exit region. The test section walls 
are built of Individual insulated copper segments to provide a uniform surface 
temperature on the segment, and each segment is provided with an individual 
foil -heater on its rear-face. Figure 4.4 shows the segments and the five 
regions of uniform temperature that result. The top and bottom walls are 
mirror-images of each other. 


The apparatus was allowed to reach equilibrium before readings were taken. Air 
mass flow rate, air pressure drops, electrical power Inputs, and segment 
temperatures were recorded. Heat losses to the background were determined for 
each configuration tested. 
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Figure 4.4 Definition of Heated Test Sections for Heat Transfer Measurements 




5.0 COMPARISONS OF CALCULATIONS AND MEASUREMENTS 


5.1 Effects of Rotation 


5.1.1 Computational Set-Up 

The computer storage capacity available at the time the calculations were made 
was represented by a 4 Megabyte virtual machine. For this machine, the code 
was configured with a 30x20x20 maximum array size. Almost the maximum allowable 
grid was used within this capacity, 30x20x14. Cells in the x and z-coordinate 
direction were arranged to have an aspect ratio of unity; in the y-coordinate 
direction the grid lines were arranged to enhance the view obtained close to 
the suction and pressure sides of the duct, where the boundary layers of 
interest were. With this grid the solution cannot be considered grid 
independent and must be viewed with a remembrance of the numerical diffusion 
that is present. The resolution of the boundary layers cannot be considered 
adequate, and it should be recalled that the wall -function approach is being 
used. 

The measured inlet profiles at an RPM of 165 for mean and fluctuating 
velocities were used as starting conditions. The usual assumptions concerning 
hydraulic diameter were made for a length scale of turbulence. 

5.1.2 Flow Visualization 

Although the experiment did not provide any visualization of the flow, the 
code does provide such a convenience. Since a sound appreciation of the basic 
flow patterns is essential for a comprehensive understanding of flows as com- 
plex as those under consideration, this aspect of the experiment is visited 
here. 

Once again, the flow visualization technique used in the calculations is the 
random streakline. Figure 5.1 shows the calculated streaklines in a cross- 
section of the duct at two axial locations along it, equidistant from the axis 
of rotation. The formation of the counter-rotating vortex pair due to Coriolis 
forces can be observed. The strength of these vortices increases with distance 
down the duct because the sign of the Coriolis force remains the same on either 
side of the axis of rotation. Rotation of the duct causes a progressive shift 
of the vortex centers from the mid-line of the duct toward the pressure-side 
of the duct. This does not change as the axis of rotation is crossed. However, 
since the developing vortices are still relatively weak on the inflow side 
(flow toward the axis of rotation) of the duct, this shift only becomes evident 
on the outflow side (flow away from the axis of rotation). 

The results observed are in agreement with the flow visualization of Reference 
20 and the calculations of Reference 21. 

The calculations agree quite well with the measurements of mean axial velocity 
profile and indicate the almost constant boundary layer thickness on the wall 
that is the pressure side on the outflow side of the duct, and the thickening 
on the opposite wall. Since the turbulence model does not recognize the extra 
strains associated with rotation, the calculated effect is due entirely to 
secondary flow development. 
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Figure 5.1 Calculated Secondary Flow Patterns Due to Coriolis Forces With 
Rotation (Moon's Experiment) 


5.1.3 Mean Velocity Development 

A comparison of the calculated and measured development of mean axial velocity 
profiles at the duct centerline is contained in Figure 5.2. 

On the pressure side of the duct, the Coriolis forces increase turbulent 
momentum exchange in the boundary layer and cause a thinner boundary layer. On 
the suction side they inhibit momentum exchange thereby causing a thicker 
boundary layer. The secondary flow development shown in the cross-sectional 
flow visualization, also causes low momentum fluid to be transferred to the 
suction side resulting in a thicker boundary layer. This profile development 

is shown in the measurements, as Figure 5.2 shows. 

# 

The calculated velocity profiles shown in Figure 5.2 were obtained using the 
"law of the wall" approach for the boundary layers, as was described in 
Appendix A2. To ascertain exactly how well this computational convenience 
really works, it is necessary to examine the boundary layer parameters. The 
boundary layer thicknesses, defined as the distance normal to the wall to 
where the local velocity reaches 99 percent of the "free-stream 1 value, are 
compared in Figure 5.3. These thicknesses reflect the profile development 
across the duct in a quantitative manner. The agreements of the calculations 
with the measurements are good. Figure 5.4 compares momentum thicknesses, a, 
for the suction side. Again, the agreement is good. 
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Figure 5.2 Calculated and Measured Mean Axial Velocity Profiles in Moon's 
Experiment 


Figure 5.5 compares measured and calculated skin friction coefficients along 
the duct. The measurements are based on a Preston tube, and the calculations 
were obtained on the basis of a Clauser plot. It can be seen that on the suc- 
tion side of the duct there is excellent agreement between measurement and 
calculation, except for the last station where the measured skin friction is 
higher than calculated. The last measurement station is close to the duct out- 
let, and flow at this station was probably influenced by the outlet (Ref. 17). 
On the pressure side, the initial agreements are good, but the measured skin 
friction increases with distance much more rapidly than calculated. This is 
most likely due to the turbulence changes due to rotation, which are not 
accounted for in the calculation. The turbulence changes and the secondary 
flows may also causes changes in the velocity profile through the boundary 
layer. This is assumed to be "universal" in the calculations. 
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5.1.4 Fluctuating Quantities 

The turbulence model used in the calculations is known to have both limitations 
and shortcomings for the present application. Therefore, good agreements be- 
tween measurements and calculations for turbulence quantities should not be 
expected. For this reason, extensive comparisons of these quantities were not 
attempted. 

Inherent in the turbulence model is the assumption of turbulence isotropy, 
i.e., if a prime represents a fluctuating velocity component, u 1 =v 1 =w ' . The 
code calculates the specific kinetic energy of turbulence K (Equation 3.3) 
which, with isotropy, is: 


K = 3/2 (u ' 2 ) 


(5.1) 


where the overbar denotes a time-mean velocity. Now, the experiment presents 
measurements of u 1 and v'. If an assumption is made for w 1 , the "measured" and 
calculated specific kinetic energies of turbulence may be compared. There are 
two plausible assumptions for w 1 : 

1. w'=v', so that K = 1 /2 ( u ' 2 + 2 v' 2 ) (5.2) 

2. w'=1.5v', so that K = l/2(u' 2 + 2.5 v' 2 ) (5.3) 

The latter assumption is based on the flat plate boundary layer measurements 
of Klebanoff (Ref. 22). 
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Figure 5.6 compares the calculated profile of K across the duct at the 66 inch 
n fift ml station with "measured" profiles based on the measured profiles of u 
lift lo etJer wU* ttTtwo assumptions above forW. ““tside of the boundary 
lavprs the three curves are in good agreement. The agreements with either of 
Iff? assumptions U are S quite poor in the thicker boundary layer existing on 
the suction wall. The agreements were better for the thinner. • 

boundary layer, especially for the second of the assumptions about w . 
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Figure 5.6 Kinetic Energy of Turbulence Profiles Across Moon's Duct Close 
to Exit 


The Reynolds stress u'v 1 can be found (Equation 3.2) from: 


-u‘ v' 


* au + 

- ay ax' 


(5.4) 


The calculated mean velocity information and the turbulence model (Equation 
3.4) provide the gradient information and eddy viscosity to enable the Reynolds 
stresses to be found. These are compared with the hot-wire measurements of the 
Reynolds stress u’v 1 in Figure 5.7 at two stations down the duct. 
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On the pressure side of the duct - u'v 1 is negative, and on the suction side 
it is positive, and this is reproduced in the calculations. However, the 
magnitudes of the Reynolds stresses are over-estimated on both sides of the 
duct by the turbulence model. The discrepancy between the measurements and 
calculations increases with increasing distance down the duct. It should be 
noted that the accuracy of the gradients used in Equation 5.4 is severely and 
adversely affected by the limited grid used in the calculations. 

5.1.5 Surma ry 

The experiment is a nice one because it confines itself to a simple geometry 
and isothermal flow, and has mean and fluctuating velocity information. 
Boundary layer parameter information is also included. It lacks flow visuali- 
zation and cross-stream velocity component measurements, and the instrumenta- 
tion is Intrusive. Measurement accuracy information Is not available. 

The calculations reproduce the measured flow behavior in every respect in 
qualitative fashion. The accuracy is affected by the relatively coarse grid 
that had to be used. The quantitative accuracy of the mean flow calculation 
is, however, fairly good. The boundary layer parameters do not calculate 
particularly accurately. Considering the grid used, the use of wall functions 
to describe the boundary layers, and the known limitations of the existing 
turbulence model, this should not be surprising. The fluctuating quantities 
are not accurately calculated. This is common experience with the two-equation 
turbulence model (Ref. 3) when used for complex flows. 

5.2 Effects of Duct Aspect Ratio with Rotation 

5.2.1 Computational Set-Up 

The available storage capacity on the computer was again restricted to 
4 Megabytes. However, since the previous calculations had required large 
amounts of computer time to achieve convergence, some economies were 
Introduced to reduce this. This of course, represents a conscious trade-off of 
resolution and accuracy for economy. The grid used was 12 x 18 x 12. The 
y-grid lines were not uniformly distributed, but expansion and contraction 
factors were applied to the grid to concentrate cells on the suction and 
pressure side walls where the boundary layers of interest are. 

5.2.2 Flow Visualization 

Figure 5.8 uses streaklines to show vortex development in the cross-sections 
of two ducts with aspect ratios 1/2:1 and 7 1/3:1, rotating at a steady speed 
of 175 RPM. The cross-sections are 1.37 m (4.5 ft.) from the inlet. 

The general features of the flow are similar to those shown in Figure 5.1, and 
have the counter-rotating vortex pair formation due to Coriolis forces 
(Figure 5.8). The effects of the passage aspect ratio, which interchanges 
major and minor axes, are obvious. 
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Streaklines Showing Vortex Development for Different Aspect 
Ratio Ducts in Moore's Experiment 


5.2.3 Velocity Profile Development 

To examine the effect of passage aspect ratio centerline profiles of mean 
axial velocity are compared at a station of 1.37 m (4.5 ft.) from the entrance 
to the duct. This is on the outflow leg of the duct. This is done for aspect 
ratios of 1:1, 4:1 and 7 1/3:1; the rotational speed was again 175 RPM. 

Figure 5.9 compares the calculated profiles. For reference, the aspect ratio 
was 2:1 for the experiment described under Section 5.1. To obtain the aspect 
ratio variation, the duct width D was kept constant and the height L was 
varied. As the aspect ratio, L/D, increases past 1:1, the profile shift to the 
pressure wall due to rotation that was described in Section 5.13, becomes 
reduced. The flow visualization of Figure 5.8 suggests this is because the 
vortex centers become further apart, and the flow on the horizontal centerline 
is not well -organized, remaining relatively undisturbed. 

The calculations show the effect of aspect ratio with qualitative accuracy. 

The quantitative agreement is reasonable for the 4:1 aspect ratio. It is not 
as good for the 1:1 aspect ratio. The accuracy for both aspect ratios is 
inferior to that for the 2:1 aspect ratio passage calculated earlier. The dis- 
agreements are particularly severe on the suction wall where the boundary layer 
thickening takes place. However, careful consideration of Figure 5.10 yields 
two relevant observations: 
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CENTERLINE VELOCITY PROFILES 
PREDICTIONS FOR J. MOORE EXPERIMENT FOR ROTATING DUCT 
FOR ASPECT RATIOS OF 1:1, 4: 1,7 1/3:1 
FOR A ROTATION NUMBER OF 0.0207(175 RPM) 
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Figure 5.9 Calculated Effects of Duct Aspect Repeat Ratio on Mean Axial 
Velocity Profiles in Moore's Experiment 


Figure 5.10 compares the measured velocity profiles with those calculated for 
aspect ratios of 1:1 and 4:1. This comparison provides for a quantitative 
assessment of the code. 

1. The normalizing factor u max was poorly defined in the experiment, and 
quite clearly the value used in Figure 5.10 is not the correct one for 
any of the profiles shown. 

2. The gradients of velocity are fairly-well calculated. 

Therefore, there is absolutely no reason to suppose that the quantitative 
accuracy is any worse, or better, than that demonstrated in Figure 5.2 for the 
2:1 aspect ratio passage. The effects of the grid changes, or the adequacy of 
either grid, cannot be assessed. 
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(RPM=175) 



Figure 5.10 Comparison of Calculated and Measured Aspect Ratio Effects on 
Mean Axial Velocity Profiles 


5.2.4 Effects of Rotation 

For a fixed aspect ratio of 1/2:1 centerline profiles of mean axial velocity 
were calculated for a number of rotational speeds covering rotation numbers 
S2 D/{j , from 0 to 0.082. The results are shown in Figure 5.11 (See Figure 
5.8 for the general flow pattern at this aspect ratio). 

As might be expected, the profile shift toward the pressure side increases 
with increasing rotation. This is in qualitative agreement with the findings 
in Reference 18; direct comparisons against the data were not made because of 
the uncertainty concerning the definition of the normalizing factor u max in 
the reference. 
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Figure 5.11 Calculated Effect of Rotation Number on Mean Axial Velocity 
Profiles in Moore's Experiment 


5.2.5 Summary 

It is unfortunate that the full potential of the experiment could not be used 
due to inadequate definitions in the data presentation. However, the calcula- 
tion, once again, reproduced in qualitative fashion the measured behavior in 
every respect for the limited parameters compared. If more comprehensive com- 
parisons could have been made there is no reason to be expect results any 
different than found under Section 5.1. The accuracy would still be dominated 
by the coarseness of the grid upon which the calculations were performed. The 
limitations imposed through the use of wall functions and the two-equation 
(isotropic) turbulence model would remain. 

The combined results show that passage aspect ratio introduces another impor- 
tant effect that influences the flow, in addition to rotation and the presence 
of a sharp bend. It might be anticipated that combinations of parameters could 
be such that ranges of variables would exist where one or another of the 
effects dominated the overall flow field. 
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5.3 Surface Heat Transfer in a Sharp Bend 

5.3.1 Computational Set-Up 

The flow path shown in Figures 4.3 and 4.4 was provided in the experiment with 
a plenum chamber at entrance and an exit section. To provide better definition 
of the inlet conditions to the section, and to give the correct back-pressure, 
both the inlet and outlet sections were calculated. A 12 x 6 x 8 grid was used 
for each of these. For the heated test section, a 24 x 14 x 8 grid was used, 
arranged to increase the number of grid lines adjacent to the surfaces of 
interest. This grid-densing did not, of course, approach that necessary to 
truly calculate the boundary layers. 

5.3.2 Comparison Conditions 

From the collected test data, an inlet width W1 (see Figure 4.3) of 3.175 cm 
(1.25 in.) was selected to give a ratio W1/W2 of unity. The duct height D was 
1.27 cm (0.5 in.) and the turn dimension, H, was 2.54 cm (1 in.). Two Reynolds 
numbers, 10,000 and 60,000, were chosen. 

5.3.3 Flow Visualization 

Figure 5.12 shows a flow visualization made in a similar geometry to that 
under consideration (Ref. 23). It serves to reveal the basic features of the 
flow in such a duct. The flow shown is the surface flow made visible using the 
technique of Langston at the University of Connecticut. 
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Figure 5.12 Experimental Visualization of Flow Development on Bottom Surface 
of a Sharp Turn, Showing Separated Flow Regions 
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The calculated surface flow (z = 0.001 ft or 0.03 cm from the surface) is 
shown in Figure 5.13. The corner-vortex on the first bend is reproduced in the 
calculations, as is the separated flow region on the dividing partition imme- 
diately downstream of the second bend. The presence of a thick boundary layer 
on the partition downstream of the separated flow reattachment appears to be 
evident. 
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Figure 5.13 Calculated Streakline Flow Visualization on Bottom Surface of 
Metzger's Experiment Showing Separated Flow Regions 


5.3.4 Surface Heat Transfer 

The Nusselt numbers calculated for the lower surface of the duct are compared 
with values derived from measurement in Figure 5.14 for the two Reynolds num- 
bers, which are based on hydraulic diameter. The Nusselt numbers are overall 
values for each of the regions shown in Figure 4.4. The calculated heat trans- 
fer coefficients were based on a "film temperature" which was taken as the 
value at the first grid-node away from the wall, a specified wall temperature, 
and the calculated heat flux. 

The results of the comparison are encouraging. The level -change in Nusselt 
number with Reynolds number is reproduced, as is the large increase caused by 
the 180° sharp bend — a doubling of Nusselt number. On the outflow leg of the 
passage, calculated Nusselt numbers are less than measured, indicating a 
faster recovery of the flow from the effects of the bend. This may be due to 
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an incorrect calculation of the separated flow region on the downstream side 
of the partition. The numerical diffusion Introduced due to the coarse grid 
could be responsible for this. A further cause could be the failure of the 
turbulence model to account for the changes in turbulence structure associated 
with the sharp bend. There is some question also as to how representative the 
temperature gradient to the wall is, and used in calculating the heat flux, 
with such a coarse grid. 


(Lower surface) 
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Figure 5.14 Comparison of Measured and Calculated Nusselt Number at Two 
Reynolds Numbers in Metzger's Sharp Turn Experiment 


5.3.5 Summary 

The surface heat transfer should be the most sensitive parameter to the short- 
comings known to exist in the current code; specifically, the lack of resolu- 
tion and the numerical diffusion associated with the forced use of coarse 
grids, and the isotropic turbulence model. There is reason to believe these 
shortcomings did adversely influence the calculations. However, the results 
are still extremely encouraging and do reproduce the general features, the 
changes with Reynolds number, and the correct levels of heat transfer. 



6.0 DISCUSSION OF RESULTS 


6.1 General Observations 

A program of verification testing such as that described above in Sections 3.5 
and 5.0 is a necessary step in the application of any CFD computer code to a 
new area; in this case, to heat transfer with rotation, which is far removed 
from the original application of the Pratt & Whitney-TEACH code to chemically- 
reacting combustor flows. 

It was first necessary to check the soundness of the modifications introduced 
into the momentum equations to account for Coriolis and centrifugal forces, 
and that they were correctly programmed into the code. It was also desirable 
to obtain a "feel" for the influence of the known limitations of the code on 
the solutions obtained. This was necessary to ascertain what might have to be 
changed later in order to improve the quantitative accuracy. Furthermore, it 
was important to build up user experience of using the code for three-dimen- 
sional heat transfer calculations and rotating flows which were areas new to 
the TEACH group at Pratt & Whitney. This experience was required to build 
confidence in the modified code and its three-dimensional post-processor, to 
acquire problem set-up experience, and so establish an accuracy-base against 
which to compare the subsequent calculations of the experiments to be performed 
at United Technologies Research Center under this contract. The verification 
study accomplished these goals. 

The known shortcomings existing in the modified code were concerned with 
numerical diffusion and flow resolution associated with the sparse grids, the 
stair-step representation of curvilinear geometries which is a consequence of 
the finite difference approach, use of an isotropic turbulence model not able 
to account for the changes in turbulence structure due to rotation and curva- 
ture, and the computationally-convenient use of "universal" wall -functions to 
represent boundary layers. Of these, only the second was not present in the 
verification testing carried out. The one believed to be most limiting in 
achieving quantitative accuracy was the first-numerical diffusion. This was, 
in the main, due to the unsuitability of the existing computer configuration 
for CFD work. Being concerned with surface heat transfer, it was also reason- 
able to expect that the turbulence model and use of wall functions might also 
exert a secondary effect on accuracy through their influence on the boundary 
layers. 

The verification testing carried out was limited to simple configurations; 
specifically, straight passages with rotation, sharp turns without rotation, 
and smooth walls in all cases. The benchmark experiments themselves were 
considered to be adequate by the standards established for such a purpose 
(Section 4.1). They enabled the major features felt to be important in the 
blade cooling passage configuration to be explored. These were rotation, 
passage aspect ratio, and sharp turns. A more extensive program of verification 
testing against a broader base of benchmark experiments would have been desir- 
able, but this desire was constrained by considerations of time and expense. 
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The general results of the verification study were exceedingly encouraging 
given the state-of-the-art in turbulent, recirculating flow CFD. The changes 
made to the equations of motion appear to be satisfactory and they seem to be 
correctly programmed. Qualitative agreements of flow patterns and data trends 
were obtained in all cases. Quantitative agreement was better than anticipated. 
The defects encountered appear to be those that were expected. 

6.2 Fluid Dynamics 

The fluid dynamic aspects of interest cover the general flow field, velocity 
profile development including boundary layer parameters, and the turbulence 
structure. 

The results of the flow visualization are satisfactory. The flow characteris- 
tics are well shown In all cases, and the pictures presented are a great help 
in understanding the developing features. The flow visualization also helps 
explain the form of some of the measurements. For example, the sharp Increase 
in heat transfer shown in the curves of Figure 5.14 can be related to the flow 
separation off the edge of the partition forming the sharp turn, as revealed 
in Figure 5.13. Similarly, the negligible shift in centerline axial velocity 
profile with rotation for the 7 1/3:1 aspect ratio passage visible in Figure 
5.9 is explained by Figure 5.8 which shows a large, relatively axial flow 
region existing between the widely-separated vortex pair. The flow visualiza- 
tion suggests that combinations of passage geometry and operating conditions 
might exist under rotation where one flow feature or another might dominate 
the overall characteristics, or the characteristics in one section of a multi- 
pass passage (e.g., see Figure 3.8). 

The development of axial -velocity profiles under the Influence of rotation is 
calculated with acceptable accuracy, showing that at least the right equations 
of motion are being solved correctly. Indeed, there would be something 
seriously wrong with the modified code if centerline profiles of axial velocity 
could not be calculated. 

The major discrepancies evident in Figure 5.10 are believed to be due entirely 
to the uncertainty with the values used for Um ax . Obviously, by definition 
all the curves must have a value somewhere of unity for u/u max . This being 
so, the effects of passage aspect ratio (on centerline profiles) may be con- 
sidered also correctly calculated. 

The code does go astray when it comes to calculating the boundary layer para- 
meters with rotation. Although the boundary layer thickness u is adequately 
calculated on both suction and pressure sides with rotation (Figure 5.3) and 
the momentum thickness a on the suction side is reasonably well calculated 
(Figure 5.4) the skin friction is not so well calculated. Figure 5.5 shows 
good agreement on the suction side, but serious disagreement on the pressure 
side. The implication of this is that the universal profile used in the wall 
functions (see Appendix A2) is a good approximation on the suction side of the 
duct, but that the actual profile on the pressure side is different. 
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Moon (Ref. 17) reported that his measured boundary layer profiles were differ- 
ent on the two sides of the duct; the profile on the suction side was better 
represented by a power law. Moore (Ref. 18) showed that the measured departures 
from a log-profile on the suction side are functions of both rotation number 
D/“, and passage aspect ratio, with the largest effects being for aspect 
ratios of unity and less. 

The effects reported by Moon and by Moore are related to the suction side of 
the duct, while the discrepancies in skin friction revealed in Figure 5.5 are 
related to the pressure side of the duct. The measured skin friction was 
obtained from Preston tubes and was reported as being about 10 percent higher 
than that obtained from Clauser plots (Ref. 24). The 10 percent difference 
between the Clauser plot results and the Preston tube measurements is far less 
than the discrepancy between the calculated and measured skin friction on the 
pressure side. At the beginning of the duct, the measured and calculated skin 
frictions are in agreement, so it does appear that the calculation is failing 
to account for all the physics associated with the developing flow in the duct 
under the influence of rotation. 

A comparison of measured and calculated centerline profiles for Moon's 
experiment is made in Figure 6.1 at 106.7 cm (42 in.) from the duct origin. 

The information is presented in Clauser plot form so that the relationship to 
skin friction may be readily appreciated. The universal relationship used is 
stated on the figure; for reference, the TEACH code uses a relationship that 
can be expressed as: 



from which the calculated skin friction values given in Figure 5.5 were 
obtained. Equation (6.1) yields slightly lower values of skin .friction than 
the equation in Figure 6.1. 
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UNIVERSAL PLOT OF VELOCITY PROFILE AT STATION 4M = 165 RPM 
Figure 6.1 Clauser Plot of Axial Velocity profiles for Moon's Experiment 


Although the absolute value of skin friction from Figure 6.1 cannot be com- 
pared directly with those in Figure 5.5 since the calculations used a different 
universal profile. Equation (6.1), relative magnitudes may be compared. 

It can be seen from Figure 6.1 that although calculated and measured profiles 
for both suction and pressure sides of the duct agree with each other in the 
outer "wake" component of the boundary layer, this is not so in the log-law 
region of the inner layer. On the pressure side, the calculated profile indi- 
cates a lower skin friction at this station than the measured profile, while 
on the suction side the calculated profile indicates a higher skin friction 
than measured. Both calculated and measured profiles indicate that skin fric- 
tion is higher on the pressure side than on the suction side, and this is 
consistent with the Preston tube measurements shown in Figure 5.5. The differ- 
ences between calculation and measurement should not be as large as that 
indicated for the pressure side in Figure 5.5. Comparison of the calculated 
suction side profile with the measured profile shows that the calculated pro- 
file is indicating the logarithmic form which gave rise to it, while the 
measured profile is not logarithmic for the smallest value of y U v of about 
2,500. 
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The Clauser plot approach to skin friction implies that a universal profile 
has a fixed turbulence structure associated with it, i.e., 


C 



T w / U 2 


V 


( 6 . 2 ) 


and. 


t w = -pu'v' (6.3) 

The skin friction calculations are therefore wrong in two respects: the wrong 

velocity profile is calculated because of the imposed logarithmic wall func- 
tion, and the isotropic turbulence model does not calculate the changes in 
turbulence structure associated with rotation, specifically the increase in 
turbulence on the pressure side and its suppression on the suction side (see 
Figures 5.6 and 5.7). 

6.3 Heat Transfer 

The agreement of calculated Nusselt numbers with those derived from measure- 
ments at both Reynolds numbers investigated are satisfactory in the inflow 
passage and around the sharp turn. They go astray in the outflow passage where 
the calculation shows lower heat transfer than measured. However, the form of 
the heat transfer distributions and the increases with Reynolds number, are 
correctly calculated. Figure 5.14. 

The implication of Figure 5.14 is that the calculations are showing a much 
faster recovery from the effects of the bend and a return to passage flow than 
is measured. This is almost certainly associated with calculation of the 
separated- flow region on the downstream side of the central partition forming 
the turn. This separated flow is shown in Figures 5.12 and 5.13. The major 
cause of the failure to calculate the separated region correctly is almost 
certainly a lack of resolution and the high numerical diffusion arising from 
the coarse grid used. In addition, the "law of the wall" used in the wall 
functions contained in the code is unlikely to hold in the vicinity of the 
reattachment point (see Appendix A2 for details). 


6.4 Summary 

Although results of the study indicate an encouraging start in calculating 
heat transfer in rotating passages involving sharp flow-turns, they also in- 
dicate areas where the calculation is deficient. The difficulties revealed 
were not unexpected (see Sections 3.2 and 3.3). Specifically, it is clear that 
the two-equation turbulence model is not adequate, and representation of the 
extra strains associated with flow rotation and turning is needed. The use of 
wall -functions might be acceptable if these were modified to account for rota- 
tion. Resolution and numerical diffusion with the allowable grids (computer 
and cost limitations) prevents accurate determination of the gradients near 
the walls. The lack of a curvilinear formulation for the equations was not a 
factor in the present study, but could create a problem in the experiment to 
be calculated (Figure 3.3). Ironically, the calculation is not so likely to be 
limited by this lack for real blade passages (Figures 2.1 and 3.5). 
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7.0 CONCLUSIONS AND RECOMMENDATIONS 


7.1 Conclusions 

As a result of analysis and data correlation in Phase I, Task III of the NASA 
sponsored program to investigate heat transfer within rotating turbine blades, 
the following conclusions may be drawn: 

1. The 3D-TEACH code is a suitable selection as a baseline code. 

2. The modifications to the momentum and enthalpy equations to account 
for rotation are appropriate. 

3. The accuracy of the code for making three-dimensional calculations 
where the "mainstream" and "boundary layer" regions of the flow are 
both impor- tant, is severely compromised by the relatively coarse 
grids that have to be used with the current generation of computers. 

4. The turbulence model currently incorporated in the code does not 
permit adequate calculation of the turbulence quantities where 
rotation and turn- ing is involved. 

5. The lack of a body- fitted coordinate system in three-dimensions causes 
difficulties in adequately representing smooth turns. 

6. Where a stair-step representation of the passage walls is important, 
the existing wall -functions could be Inappropriate. 

7. The results of the study are encouraging enough both to believe that 
CFD can make a useful contribution to blade internal cooling, and to 
proceed further with the present program. 

7.2 Recommendations 

Following-on from the conclusion, it is possible to make some recommendations. 
It is not implied that these recommendations will be Implemented during the 
course of the present contract. 

1. Turbulence models more advanced than the present two-equation model 
can be formulated and should be explored for this application. 

2. Modification of the wall functions to account for rotation can be made. 

The computer problem and its resolution is beyond the scope of this report. 
Development of a three-dimensional curvilinear code is a difficult and costly 
task. It is certainly beyond the terms of the present contract. 


48 



APPENDIX A1 


DESCRIPTION OF 3D -TEACH 
1.0 Description of Aerothermal Model 

3D-TEACH is a computer code that can solve fully three-dimensional fluid 
dynamics problems. It can handle axisymmetric, planar, or three-dimensional, 
elliptic, turbulent flows. It is one of a family of such codes with titles 
2D-TEACH, 2D(C) -TEACH, 3D -TEACH and, 2D-PREACH. 

The input to all of these codes is generalized such that different problems 
can be run without the need for Fortran programming between problems. Also, 
the physical models used can be turned on or off by input command. These 
collective features result in an extremely flexible system. 

The codes form a family in that: 

a) As computer codes they are written with the same format, menus and 
commands such that an operator trained to run 2D-TEACH can easily run 
3D-TEACH. 

b) All codes have an interactive nature using the IBM Conversational 

Monitor System (CMS), and use prompts and cautions to ensure smooth 

execution of a case. The operator has the choice of either CMS or 
Batch running. The recommended procedure is to set up a case and get 
it running on CMS, then switch it to Batch to complete execution. 

File modification and selection of running mode is identical for all 
codes. 

c) The same basic equations are solved and the same physical models are 

used, together with solution algorithms. In all codes such that 

regression is possible. This means that the same two-dimensional 
problem can be solved on 2D-TEACH and 3D-TEACH, and the same results 
will be achieved. The only difference between the 2D-TEACH and 
3D-TEACH is the additional dimension available In the 
three-dimensional code. Also, the post-processors available with 
3D-TEACH are necessarily more comprehensive than those with 2D-TEACH, 
in order to adequately display the results. 

The acronym TEACH (Teaching £lliptic Axisymmetric Characteristics 
Heurlstically) represents a generic solution technTque and these codes 
represent current production state-of-the-art calculations In terms of 
equations solved, physical models used, discretization of the equations, and 
solution algorithms. They are not perfect, but are a marked advance on 
one-dimensional flow calculations and global modeling that formed the previous 
capability. The structure of the codes has been made such that modular 
replacement can be carried out as better models and solution algorithms are 
developed, while the basic framework and operational features remain. 

The 2D-PREACH code is similar in concept to the others, but uses different 
solution procedures. It is not considered further here. 
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1.2 Outline of Calculation Procedure 
1.2.1 Equations to Be Solved 

The combustion process is an extremely complex turbulent flow. It is a 
somewhat daunting task to set about describing such a random flow of 
chemically active eddy structures in terms that can be currently solved and 
can provide useful answers for the designer of practical equipment. 

Currently, the most practical approach is to stay within the framework of 
continuum mechanics and to use a statistical description of the turbulence, 
coupled with the accepted Eulerian description provided by the Navier-Stokes 
equations of motion. Hence, an Instantaneous quantity is described as the sum 
of a time-averaged value and a random, fluctuating value. 

When the statistical description of an Instantaneous quantity is substituted 
into the Navier-Stokes equations (Ref. A1 ) and time averaged, the resulting 
equation set Is known as the Reynolds equations (Ref. A2). These equations are 
similar to the Navier-Stokes equations except that time-averaged quantities 
are used, and for the appearance of time-averaged correlations of fluctuating 
quantities. 3 

Turbulent motions increase the apparent viscosity of a fluid by some orders of 
magnitude since there Is a continuous transfer of energy from the mean flow 
into large eddies and thence, cascading down through progressively smaller 
eddies, to the molecular level where the energy Is dissipated as heat. If 
laminar diffusion terms are therefore very much smaller than turbulent 
diffusion terms, then neglect of fluctuations in laminar viscosity is 

permissible. Th* 

frequently used 
density, althou' 
not large. This 
equations. 

The simplified Reynolds equations are expressed in terms of time-mean 
quantities and also cross-correlations of fluctuating velocities such as 
PUf Uj'. These terms are known as the Reynolds stresses, and result in a 
closure problem. Turbulence modeling provides the necessary descriptions of 
the Reynolds stresses in known or determinable quantities. When the flow 
consists of more than_o ne che mical species, modeling Is required also for the 
turbulent mass flux pu^'m-j . These terms arise from applying the 
statistical treatment ot turbulence to an Instantaneous species transport 
equation. The Instantaneous energy equation is given the same treatment. 

To model the turbulent mass fluxes it was assumed that, similar to molecular 
Schmidt and Prandtl numbers, there are turbulent Schmidt and Prandtl numbers 
that relate turbulent mass and heat diffusivities to momentum diffusivity. 
Closure to the Reynolds equations was provided by a particular turbulence 
model known as the two-equation or K-e model. It relies on the eddy viscosity 
concept. Finally, the modeled equations were algebraically manipulated into a 
general form in cylindrical co-ordinates: 


is result m simp n ricaxion or tne Keynoids equations. It is a 
practice (Ref. A3) to also neglect terms involving fluctuating 
gh this implies that temperature differences in the flow are 
practice also results in simplification of the Reynolds 
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(Al.l ) 


{L (pu 0) ♦ 9 (rpv •) ♦ 9 (p W 9 ) - £ 

r^r r5« £x 





r eff,f 


±h • s# 

#0 


where: 

0 = any of the independent variables 

ef f 0 = an appropriate turbulent exchange coefficient, depending on 
what i represents 

S 0 = a so-called "source term" which lumps together all other 
terms in a given equation not included in the first four 
terms of Equation (Al.l). 

The equation given is for steady state flow. Reynolds averaging does result in 
the equations retaining time-dependent terms, but these havp been dropped. 

This was done for two reasons: (1) compatibility with the present design 
system and (2) time averaging precludes dynamic behavior other than that 
induced deliberately through one of the independent variables. 

By way of example. Figure Al-1 gives the values of some of the items in 
Equation (Al.l). 

1.2.2 Numerical Approach to Equation Solution 

The simultaneous set of main and auxiliary equations to be solved in a 
turbulent reacting flow with a liquid-fuel spray contains a significant number 
of individual equations, most of which are either ordinary or partial 
differential equations, and which are nonlinear. Numerical solution of these 
equations is necessary. Rearrangement into the general form represented by 
Equations (Al.l) and (A1.2) enables one solution algorithm to be used for all 
equations. 

Conventional numerical methods available to solve equations of these types can 
be broadly divided into finite difference and finite element methods, although 
the dividing line is not distinct. Finite differences have a considerable 
background, and most solution approaches utilize this method. 

The finite difference analog of the differential equations is obtained by 
overlaying a computational mesh on the flow domain, and obtaining the basic 
finite difference form of the partial derivatives for every r.ode of the mesh 
from a control volume approach (Ref. A4). The finite difference expressions, 
when substituted back into the differential equation?, yield a set of 
linearized, algebraic equations for every node of the mesh. Thus, there are as 
many sets of equations as there are nodes in the calculation domain. These 
sets, along with the problem boundary conditions, can then be solved to give 
solutions for the entire flow field. 
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Figure Al-1. Summary of Some of the Equations Solved In 3D-TEACH 
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ORIGINAL PAGE IS 
OF POOR QUALITY 

Standard numerical techniques can be employed to solve the finite difference 
forms of the differential equations (Ref. A5). A steady-state implicit 
solution method is often used (Ref. A6); an initial guess is made of the field 
variables, and these guesses are iteratively updated until the solutions have 
converged. Convergence is deemed to have been obtained when the absolute sums 
of the residuals of each variable over the whole grid goes below a specified 
value. 

The relevant equations describing the flow motions, the physical models, and 
the solution techniques are assembled into the computer codes to carry out 
direct flow simulations on high-speed, large-core digital computers . 

Figure Al-2 presents a flow diagram describing the calculation procedure, 
showing the assembly of the equations, the utilization of physical modeling, 
the computer solution, and the output for design use. The 3D-TEACH code 
conforms to this organization. 



Figure Al-2 Flow Diagram of Calculation Procedure 


1.3 Solution Procedure 

With reference to Figure Al-2, assembly of the equations governing the problem 
has been briefly described. The details of the computer solution of the 
resulting equation sets are to be described. 

Rearrangement of the equations into the general form represented by 
Equation (Al.l) enables one solution algorithm to be used for all equations. 
The equation set is solved using a steady state, Implicit, finite difference 














numerical procedure. Initial guesses are made of the field variables, and 
these guesses are iteratively updated until the solutions have converged. 
Convergence is deemed to have been obtained when the absolute sum of the 
residuals over the whole grid of each variable goes below a specified value. 

1.3.1 Discretization of the Equations 

The finite difference analog of the difference equations is obtained by 
overlaying a computational mesh on the flow domain to be calculated, and 
obtaining the basic finite difference form of the partial derivatives for 
every node of the mesh from a control volume approach, (Ref. A4). The finite 
difference expressions, when substituted back into the differential equations, 
yield a set of linearized, algebraic equations for every node of the mesh. To 
demonstrate, first in two-dimensions. Figure Al-3 illustrates the mesh and the 
control volume established about a considered node, P. The control volume 
approach is based on the satisfaction of macroscopic physical laws such as 
conservation of mass, momentum and energy. The conservation property is 
essential when combustion is taking place. 



Figure Al-3 Control Volume for the Finite Difference Scheme 


The code is written in both cylindrical and cartesian coordinates. The grid 
system consists of a set of coordinate lines intersecting in the x-r, x-e and 
r-e planes for the cylindrical system. In the cartesian system the grid is 
formed by the intersection of x-y, x-z, and y-z plane lines. The intersections 
of these lines form the grid nodes at which all scalar properties are stored. 
Vector quantities are stored midway between nodes. Figure Al-4 gives the 
finite difference grid control volumes for the scalar quantities and storage 
locations for the velocities in both coordinate systems. Note that compared to 
Figure Al-3, there are two additional neighboring nodes, F and B, denoting 
Front and jlack nodes in the z ore direction. The faces of the scalar control 
volume are denoted by lower case letters. Figure Al-5 shows typical scalar 
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control volumes in perspective and gives the face areas and volume. Since the 
velocity components are located midway between the grid nodes, the control 
volumes for velocity components are formed by planes passing through the gird 
lines. Note that since the control volumes for the velocity components are 
staggered (Figure Al-6), the areas and volumes for these control volumes will 
be different from those of the scalar control volume. 

The finite difference form of the general partial differential equation is 
derived by supposing that each variable is enclosed in its own control volume, 
as Illustrated In Figures Al-3 to Al-6. The general i transport equation has a 
source term %i. This Is expressed in linearized form and integrated over the 
control volume. The remainder of the transport equation Is also integrated 
over the control volume, and added to the integrated source term. This yields. 


C E *e " c W^w **■ C N ^n - $ s ♦ C B $5 - Cp • Dp ( <*£- $ p ) - Dy( $ p - 

♦0 N ( ^N” ^p) • ®S ( *p - ♦ 0 B ( $ B - ^ p ) - Op ( - $f) ♦ (Sy ♦ S p d> p ) 

(A1.2) 

In the above equation the convection coefficients are defined as, 

Cg • (p u) e * e ; Cp ■ (p w) f a f etc., 

and the diffusion coefficients are defined as: 

0£ - ( ) * ** ; op • (- reff » ’ *f 

\ A x / e r a fl ' 

etc. are the areas of the cell faces 


Certain weighting factors are introduced into the variation of 0 , the 
variable being calculated, and with the help of continuity. Equation (A1.2) 
can be manipulated and normalized to give the form, 

Ap*p • An *N * As *$ ♦ A N *W + *F*f * A B*B 4 $u (A1.3) 


where, 

Ap ■ A n ♦ A S ♦ A£ * ♦ A N ♦ Ap ♦ A B - Sp 

and Equation (A1.3) is the finite difference equation for and the main 
coefficients are defined as. 

An ■ 0 when P < -2 

An ■ On - JEfil when -2 < P_ <2 
2 


an w Cn 


when 
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Similarly for A$. A^ etc., where the cell Peclet number Is defined as, 

P • C|i/Djj tic. 


There are several differencing schemes that can be used to evaluate the 
weighting factors. The values of the coefficients A^, Ac, etc., above were 
obtained using Spaldings Hybrid Differencing Scheme (HD*) , of Ref. A7. 
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CONTROL VOLUME IN THE CARTESIAN SYSTEM 



CONTROL VOLUME IN THE CYLINDRICAL SYSTEM 


AREA 

CARTESIAN 

CYLINDRICAL 

REMARKS 

•w* ■ 

(AYAZ)| 

(ArA0r|| 

Cartesian systam can 
ba obtain ad from 


(AYAZlj 

(A rA6r)^ 

cylindrical systam by 
putting r - 1. 

•"* * 

(AXA2)| 

(ArA6r n )| 

M - 62. 

M * 

(AXAZ) + 

IAxA0r a )^ 


■f* - 

(AXAY)j 

(AxA r 


•t* - 

(AXAYty 

(A*Ar)| 


Vd | - 

IAXAYAZ)* 

(AxArA0r)+ 
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The hybrid differencing scheme Is unconditionally stable and the solution Is 
bounded. It uses second order central differencing for convection and 
diffusion fluxes when the absolute value of cell Peclet number is less than or 
equal to two. When Peclet number is greater than two, first order upwind 
differencing Is used for convection fluxes, and diffusion fluxes are neglected 
altogether. The switch of differencing is done both locally and directionally 
In the computational grid. Peclet number defines the relative importance of 
convective and diffusive transport. 

The finite difference Equation (A1.3) derived in the previous section could be 
used to obtain the velocity if the pressure field were known a priori. Since 
the pressure field is unknown, an iterative solution procedure, SIMPLE 
(Ref. A8) is used. SIMPLE is an acronym for Semi JmpUcit Method for Pressure 
Linked Equations. 

The essence of SIMPLE is that a pressure field is guessed, velocities are 
calculated from their finite difference equations, then the pressure and 
velocity fields are updated using a "pressure correction" equation which 
satisfies continuity. The procedure is repeated until the momentum and the 
continuity equations are adequately and simultaneously satisfied. The pressure 
correction equation can be derived from the continuity and momentum equations; 
the procedure is described below. 

The finite difference form of the continuity equation can be written as: 

(p ii) e a e - (p u) w ♦ (P v) n « n - ( p v) s a s ♦ (P")f a f - (pw) b a b « 0 

(A1.4) 


The momentum equations can be written as: 

A P V ’ A N V 4 A S V S* 4 A E V 4 A w V A r *F* ♦ A B v B* ♦ a s f p S* - Pp*) 

Ap Up* - A m u n * ♦ A $ u s * ♦ A e u e * ♦ Ay u y * ♦ Ap up* ♦ A B u B * + a w (P H * - P p *) 

A p w p* “ A N w N* + A S W S* * a E ♦ V^w* + Apwp* + a B«B* + *b ( p F* ' p p*) 

In the above equations the pressure term has been separated from the source 
term and the (*) superscript denotes the values obtained from solving the 
momentum equations using the guessed pressure. An incorrect guess will give 
rise to a "mass source", Mp, in each cell because the continuity equation will 
not be satisfied. The mass source can be found by using Equation (A1.4). Hence 

Mp • (p u*) e a e - (p u*) w **+ (P v*) n a n - (p v*) s « s -(?w' ) f a f + (?w') b a b 

If the above equation Is subtracted from Equation (A1.4) 

-M P - ( p o' ) e » e - ('p u') w a^ ♦ (p v‘) n a n - ( P V ) s *s *^‘)f«f * (Pw*)b*b 

(A1.5) 
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where 


Ug’ » (u - u*)j etc. 


The above velocity corrections can be calculated from the linearized momentum 
equations. 


Ap up* • a* (Py’ - Pp') 


(A1.6) 


Note that Up' in the momentum equation control volume is u w ' for the 
continuity control volume and similarly for Vn', etc.. On substitution of 
Equation (A1.6) in Equation (A1.5) and simplification. 


Ap Pp* 


♦ *S P S' ♦ V ♦ *£ V * k f r f' * V B ' 


(fll.7) 


where 


Ay ■ Uy/Ap) 

Sg - - Mp 

Ap * Ay ♦ A$ ♦ Ay ♦ Ag + Af+Ag 


Equation (A1.7) is called the pressure correction equation which Is solved to 
obtain corrected velocities and pressures, 

“e " ♦ u e' 

P e " V ♦ P' 

etc. 

The difference equation for P' (Equation (A1.7)) Is In the same form as the 
difference equations for 0 (Equation (A1.3)) and hence a single solution 
algorithm can be used to solve all difference equations embodied in the 
numerical method. 

Since the SIMPLE procedure computes the variable fields successively It is 
highly flexible with respect to the methods of solution which it will admit 
for the difference equations. At present the following line by line iteration 
method Is employed. This method Is also known as Alternating Direction 
Jjnpllclt Method (Ref. A9). The ADI methods were Initially formulated for 
unsteady equations; their adaptation to steady state equations is sometimes 
also known as Alternating Direction Jteratlve Methods. 

The finite difference Equation (A1.3) to be solved Is 

Ap Pp - Ay #y ♦ A $ #$ + Ay #y ♦ A £ # E + + A B 0 B + S u 
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where 0 stands for u, v, p, K, e , and H successively. This equation can be 
recast In the following form 


Ap *p -Ay P n ♦ A S ♦ c ' 


or 


Ap Pj - Ay Pj*i ♦As Pj.i ♦Cj’ 

To solve the equations for points on each line (e.g., N-S line) values on 
neighboring lines are assumed to be temporarily known. The equation for each 
point on the N-S line then reduces to one where only three values ( 0 p , 

0N, 0 s in Equation (A1.8)) are unknown. An equation of this type can then 
be solved by the Tri-Diagonal Matrix Algorithm { TDMA) , which is explained 
bel ow. 

Equation (A1.8) can be rearranged for the j th point as 

Pj ■ Bj pj ♦ 1 ♦ Cj Pj-i ♦ Oj 


where 

Bj • a n/ a p • c j * A S/ A P 

Oj • (Ay *y ♦ Ag *£ ♦ Af$p + AgPg + 


The points on the computation grid range from 1 to Nj in the N-S direction 
with points 1 and Nj on the boundaries. Since the boundary values 0 i and 
0 Nj are known, equations for 0 2 to0yj_i are solved. The set of 
equations then becomes: 

# 2 - B 2 *3 4 c 2 *1 4 °2 

*3 " b 3 ♦u * c 3 *2 * °3 
• • • 


*Nj-l 


B Nj-l *Mj 4 C Nj-l ♦NJ-2 + °Nj-l 


(A1.9) 


Now since 0i is known 02 can be eliminated from Equation (A1.9) and so on, 
yielding a general recurrence relation 


♦j ’ A J # J*1 4 D j‘ 


(AT .10) 


To get the relation for Aj and Dj ' Equation (A1.10) is written as 

♦j-1 ■ Aj.i ♦ Oj.1 
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Now putting in the value of 0j 


■(. — -is — ) * [ fiHiLfkl 

\ A P - A S • A j-1 / j L A P - A S A j-1 J j 


(Al.ll) 


Comparing Equation (A1.10) and (Al.ll) yields coefficients for the recurrence 
f ormul a 


where 


A j " (a#/ <Ap " A S A J-1 0 J 


Dj * (< A S °j-l + c j ) / < A P * A S Aj-1 ) ) j 


(A1.12) 


(A1 .13) 


Using Equations (AT .12) and (A1.13), 0 j can be calculated from 
Equation (A1.10). Having solved for 0 j on one N-S 0 j« § on the next N-S 
line are solved and so on until the entire solution domain is swept. The same 
treatment is then applied in the W-E direction and finally, in the F-B 
direction. It is usually necessary to sweep between 1 and 3 times per 
iteration for optimunj solution time. 

The coefficient matrix formed by the finite difference equation of 0 should 
satisfy the stability condition, 

A P - „ I A n | 


Now 


A P ‘ n A n - S P 


So if 


Sp < 0 


the above criteria is satisfied. In the solution procedure care is taken so 
that Sp is always less than or equal to zero. 

In the process of the computations, convergence is assessed at the end of each 
iteration on the basis of the "Residual Source" criterion. The residual source 
R 0 is defined as 




A? 



S 
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It is required that: 


IlR^I < « **Ref 


for each finite difference equation. 

R 0 Ref is the fixed flux of the relevant extensive property fed into the 
domain of calculation, and e is of the order of 10 - 3 . 

When it is the equations for mass fraction of species that are being solved an 
additional convergence criterion requires that the sum of the mass fractions 
at each node is < to (1 + € ). 

When the flow is of variable density it is initially required that the change 
in density in one iteration at every node must also be less than e 



Since the enthalpy values in the calculation domain do not conform with the 
specie mass fractions during the first few Iterations, temperature and density 
are not updated for the first 10-25 iterations. If the density gradients are 
steep, density is updated every second or third iteration after the first 
update. 

A typical convergence plot is shown in Figure Al-7. 

Since the finite difference equations are nonlinear in nature, the convergence 
is facilitated and' sometimes divergence is avoided by under-relaxing the value 
of being calculated as: 


R . f Hew . (i-F) Old 
♦p *p # P 


(A! .14) 


where F is an under-relaxation factor which is less than one. 

The way in which the above relation is introduced into the numerical procedure 
is as follows: 


A^ » Ap/F 

s; - S u Ml - F) Af .Old 
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Figure Al-7 Typical Convergence Plot 


It can easily be shown that the effect of introducing the above modifications 
is to under-relax 0 p according to Equation (A1.4). From Equation (A1.3) we 
have 


New - n * n 2 Su 

♦p Ap (AT. 16) 


putting in the under-relaxation factors. 


* 


R 

P 


5 *n ♦ S U R 

A« 

P 


putting the value of A p R and S U R from Equations (A1.14) and (A1.15) 
in (A1.16), 


64 



I*fl ♦„ ♦ S„ ♦ (1-F) AR ,01, 

R . » P P 

*P A R 

P 

R n (n An * n * Su \ F ♦ M-F) 01d 

♦p \ — VW * p 


F . New p + (l -F) 01 d 

$p 'Pp 


It should also be noted that the effect of under-relaxation is to make the 
coefficient matrix more diagonally dominant. 

The various steps in the numerical procedure can now be summarized as follows: 

1. Guess fields for all variables. 

2. Assemble coefficients of momentum equations and solve for U* and V* 
using prevailing pressures. 

3. Solve the pressure correction equation and update velocities and 
pressures. 

4. Solve equations for other variables. 

5. Update fluid properties such as viscosity and density. 

6. Test for convergence. If not attained use prevailing fields as new 
guesses and repeat from step 2 until convergence is attained. 

In general, it is necessary to specify 0 or its gradient at the boundaries 
of the calculation domain. There are six types of boundaries: 

1 . Axi s of symmetry 

2. Unspecified wall 

3. Speci fi ed wal 1 

4. Unspecified opening 

5. Specified opening 

6. Specified blockage 

A specified boundary is one for which all boundary values such as velocities, 
temperatures, etc. are given. An unspecified boundary is one for which the 
boundary values are calculated by the code. 
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On an axis of symmetry the gradient of all 0's except v-velocity is put to 
zero; v-velocity itself is set to zero. Most walls are specified, with all 
velocities set to zero (no slip condition). A moving wall is modeled with the 
no slip condition by specifying nonzero velocities in the plane of the wall. A 
porous wall is modeled by specifying nonzero velocities normal to the plane of 
the wall, or by inputting the wall as unspecified. At the outflow, it is 
required that there be no negative axial velocity components. At high Reynolds 
numbers this requirement makes specification of boundary values of all 0's 
except u redundant. The axial velocity is specified thus, 

u "1. j " u ni-l , j ♦ Uinc 


where ni is the outflow boundary and is calculated such that the total 
mass outflow is equal to main inflow. A1 ternatively , if an exit velocity 
profile is known, it can be specified. 

The calculation mesh is constrained by the coordinate system, which presently 
has been selected as orthogonal. Therefore, curvilinear geometries have to be 
represented in the form of discrete steps or "staircases. " The specified 
blockage boundary condition permits this representation and allows inflow and 
outflow through elements of these staircases. In addition, the condition 
allows solid bodies to be placed inside the flowfield and to contain mass 
sources or sinks within them. This capability is written in generalized form 
and confers considerable geometric flexibility on the code without the need 
for Interproblem reprogramming. 

Adjacent to solid boundaries, the local Reynolds number of the flow based on 
local velocity and distance from the wall becomes very small and the 
two-equation turbulence model, which was developed for high Reynolds numbers, 
becomes inadequate. Although a version of the two equation model that can 
handle both high and very low local Reynolds numbers exists (Ref. A10), its 
application requires a large number of grid nodes (more than 30) in the wall 
layer. This is due to the steep gradients of properties in the wall region 
(Ref. All). 

To avoid these difficulties, it was argued that the flowfield in the 
calculation domain is not influenced to first order by the details of the flow 
at the walls (Ref. A12). Consequently, as a matter of computational efficiency 
and economy, the high Reynolds number version of the turbulence model was 
retained and a Couette-flow analysis was used to give an equilibrium boundary 
layer on all solid surfaces bounding the calculation domain. The resulting 
wall functions are used to link the walls to the near-wall nodes of the finite 
difference grid. 
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APPENDIX A2 


TREATMENT OF WALL BOUNDARY LAYERS 
IN BASELINE CODE 

The calculation mesh is constrained by the coordinate system, which presently 
has been selected as orthogonal. Therefore, curvilinear geometries have to be 
represented in the form of discrete steps or "staircases." The specified 
blockage boundary condition permits this representation and allows inflow and 
outflow through elements of these staircases. In addition, the condition 
allows solid bodies to be placed inside the flowfield and to contain mass 
sources or sinks within them. This capability is written in generalized form 
and confers considerable geometric flexibility on the code without the need 
for Interproblem programming. 

Adjacent to solid boundaries, the local Reynolds number of the flow based on 
local velocity and distance from the wall becomes very small and the two- 
equation turbulence model, which was developed for high Reynolds numbers, 
becomes inadequate. Although a version of the two equation model that can 
handle both high and very low local Reynolds numbers exists (Ref. A2-1), its 
application requires a large number of grid nodes (more than 30) in the wall 
layer. This is due to the steep gradients of properties In the wall region 
(Ref. A2-2). 

To avoid these difficulties, it was argued that the flowfield in the 
calculation domain is not influenced to first order by the details of the flow 
at the walls (Ref. A2-3). Consequently, as a matter of computational 
efficiency and economy, the high Reynolds number version of the turbulence 
model was retained and a Couette-flow analysis was used to give an equilibrium 
boundary layer on all solid surfaces bounding the calculation domain. The 
resulting wall functions are used to link the walls to the near-wall nodes of 
the finite difference grid. The procedure is described below for the momentum 
transfer, heat transfer, and turbulence processes. 

The wall layer is assumed to be one of constant shear stress ajj the wall 
(t= t w ). The heat flux is also assumed to be constant (q" = q w ). It 
should be noted that these conditions are true only for an impermeable wall, 
with zero or negligible streamwise pressure gradient 



The momentum equation can be reduced to a simple form as shown below. It Is 
assumed that the wall Is parallel to the x-axIs. In the case of a vertical 
wall u will be replaced by v and y by x. 

i- („ ♦ w ) lii . 0 

3y c a y 


Integrating 

(u + v c ) & l y m 0 

*y 1 0 


at y=0, fj. t =0 


so 


dy 


y“0 


(p ♦ P t ) ^ 

dy y-y 


or 


T ■ (v ♦ p ) — 
C dy 


or 


(A2.1) 


-X_ . 


(1 


where 


P dy 



Near the wall the local Reynolds number changes considerably and the approach 
adopted depends on the value of the local Reynolds number, y + , based on 
distance y from the wall and friction velocity u t. 

For convenience the wall reaion is divided Into two layers y + < 11.63, a 
fully laminar region, and y* >11.63, a fully turbulent region. Then for 
y + <11.63,Mt« . Hence from Equation (A2.1) 
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(A2.2) 


y 


An 

dy 


Also for y + > 11.63Mt>>M and t/t w = 1. 
Hence from Equation (A2.2) 

u t du + 

— — T - l 
M dy 

From log-law of the wall (Ref. A2-4), 


u + ■ 1 /k log fi (Ey + ) 


du _ 1 ' 

. ♦ _ * 
dy <y 


n t ■ My 


Hence from Equation (A2.1) 


t - t< *y* 

dy 

It should be noted that E Is an Integration constant that depends on the 
magnitude of the variation of shear stress across the layer or the roughness 
of the wall. The value of E used in the program Is for smooth walls with 
constant shear stress. Effects of mass transfer across the layer and severe 
pressure gradients can be Incorporated by modifying E which will then no 
longer be a constant. 

In the case of flow with swirl, the axial velocity u is replaced by the 
resultant velocity, ur, where 

Ug ■ + w* 

for walls parallel to the x-axis. 
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The shear stress t then becomes the resultant shear stress t r and the shear 
stress in x and ©directions can be obtained by resolving this stress 


T_ ■ 


U R 


If the wall is isothermal or the temperature distribution is specified, the 
following treatment, which is the same as the treatment used for momentum 
transfer is adopted. The energy equation can be reduced to 


where 


V - (r ♦ r t ) c p 


i!l « / r lA dT+ 

q w \U V ) dy* 


4 Pu C (T -T) 
X » Tp w 


also 


T ‘ ' [ u * * * ')] 

0 = Laminar Prandtl number 

0»t s Turbulent Prandtl number 

Now for y + < 11.63, r>>r t from Equation (A2.3) 


q" ■ r c 

P dy 


For 


r*> 11-63 r « r c 


(A2.3) 


(A2.4) 


(A2.5) 
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Therefore from Equation (A2.4) 


r. dT + 
-77- 1 
V dy 


From Equation (A2.5) 



Putting Equation (A2.6) and (A2.7) in Equation (A2.3) 

«y + C dT 
q" - — — E 

* t Ve u d T 


(A2.6) 


(A2.7) 


(A2.8) 


If the heat transfer rate to the wall q* rather than its temperature is 
specified Equation (A2.8) can be replaced by 


where Q is obtained from experiments. 
For adiabatic walls 


0 


The approach adopted for the turbulence equations Is strictly valid for the 
initial sublayer where the flow is assumed to be completely turbulent: 
y + >30, but sufficiently close to the wall so that the assumption of 
constant shear stress applies (y + < 400). In this region, the local rate of 
production of turbulence is balanced by the viscous dissipation rate . This 
local equilibrium forms the primary basis for specification of turbulent 
kinetic energy and its dissipation rate at the wall. 

In case of local equilibrium, 

“P u'v' -ili ■ pc 

, 4y 


Also, shear stress near wall (Ref. A2-5), 

t w - -p ITV - p c m 1/2 k 

(Experimental evidence suggests u'v' = 0.3 K). 
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Using the "law of the wall," 


pe 



>/\) 3/2 


where k is the constant In the log-law "law of the wall." 
Hence the wall boundary value for dissipation is 

c 3/4 3/2 

e . IU K 

<y 


The value of K at the boundary is put to zero, 
energy in the cell adjacent to the wall is 


The production of turbulence 


P 


K 


du 

dy 


iu 

dy 


The dissipation term In the kinetic energy equation is found from the averaqe 
oyer the control volume. The reason for using this technique is that is 
highly nonlinear in the vicinity of the wall. Hence, 


/ pc * 'l 


7? p c 3/4 k 3/2 


*y 


dy 


The Integral of the above expression will be logarithmic and Infinite at the 
lower limit. Hence, using the log law again, 

. JL 

dy T p Ky 


. *. p c . A p c>v« 


T • * ■/ 



la 

P 


du - p C y 3/4 K 3/2 u + 


where 


U* m y* f or y* < 11.63 


and 


u + ■ log. Ey + for y* > 11.63 

K * 
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As stated earlier, these wall functions are based on the assumption of 
constant shear stress and impermeable smooth wall with negligible pressure 
gradient. Any departure from this condition is going to cause incorrect 
prediction of wall shear stress and heat transfer. In the recirculation 
region, especially near the reattachment point, it is unlikely that the 
log-law of the wall should be applicable and the predictions in this region 
will also be incorrect. However, in spite of the above apparent weaknesses the 
two equation model does an adequate job of predicting turbulence quantities. 
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APPENDIX A3 


DERIVATION OF EXTRA TERMS GENERATED 
DUE TO SYSTEM ROTATION 


In Figure A3-1, Oq, Xg, Yg, Zg Is the Inertial frame and Ot, Xi, 

Y i. Zi is the moving frame. The motion of thejatter Is described by the 
position vector R of its origin and a vector the rotation of the moving 
frame about an axis through its origin. The components of R parallel to 
0 0 X 0 , OgYg, OgZg define the position of the moving frame and the 
components of define its rotation about the axis of Inertial frame. The 
position vector rg and rj define the position of a particle 1 P 1 in two 
frames of references. The line labelled 'path' is the track of ' P 1 seen by the 
moving observer. To him the location of its path does not change with time and 
Inertial observer sees the same line traveling through space. Therefore: 


rO = n+ * (A3. 1) 


The acceleration of P seen by the moving observer is: 


d i 2f i _ \ Vi 

dt 2 dt dt 


(A3. 2) 


and the acceleration seen by 


the inertial observer is: 


d o _ ^0 Vo 

dt 2 dt dt 


(A3. 3) 


The operator dg/dt is for differentiation in Inertial frame and di/dt is 
for the moving frame. 


The relationships between the velocity and acceleration of A seen by the two 
observers can be found bearing in mind that rn varies with time because of 
motion of P, IT varies because of the motion of Oj and rj varies because of 
the motion of A and rotation of the moving frame. This gives: 


— = — + — MBxr.) 
dt dt dt 1 


(A3. 4) 
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Figure A3-1 Derivation of extra terms generated due to system rotation 


Differencing again in the inertial frame, 

d rt d,r, d. 


d o 2? o 


V" . ^ 


dr 


dt 


dt 


d i r i 

dt 


+ — ( x r.) 
dt 1 


Vi . r1 x Vi * nx Vi 


dt dt 


dt 


dt 


dt 


(A3. 5) 


A correspondence between the operator dg/dt and dj/dt can be found from 
equations (A3.1) and (A3. 4). 

Differentiation of equation (A3.1) is the inertial frame gives: 


Vo 

dt 


Vi . V 


dt 


dt 


(A3. 6) 
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From equations (A3. 4) and (A3. 6) it can be seen that 


or 


d 0 ? 0 

V 

dj^i 

+ + (fi x r.) 

d O ? l , d 0 R 

dt 

dt 

dt 1 

dt dt 

d 0 ? l 

V r l 

+ (fiX ) 


dt 

dt 

1 



(A3. 7) 


or 


d 0 d l 

— = 

dt dt 


(A3. 8) 


Equation (A3. 8) is the relationship between the derivatives in two frames. 
From equations (A1.5) and (A1.8). 

d o 2f 

J s 




dt 


+ r 1 x 


dt 

d 0fi 


dt 



d o 2 * . d l 2f l . - d l F l . - ■‘o 

— 7 7 + fi x r i — 

dr dr dt 1 


r 

dt 


d l ? l 

+ n x + (fix r.) 

dt 1 


d o% 

“1? 


^ *5x^1. r. *2 

dt dr dt 1 dt 


- d l F l - 

+ fi x + fi x (n x ?,) 

dt 1 
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(A3. 9) 


d 0 2 R 


dt 


d i 2f i 


dt 



d l r l - - - 

x r. + 2n x + n x (six r.) 

1 dt 


The terms in equation (A3. 9) can be interpreted as follows: 


(a) 


d o% 


dt 


Acceleration of particle P to an observer in inertial frame. 


(b) 



Acceleration of moving frame (translation) as observed from 
inertial frame. 


(0 


d i 2? i 



Acceleration of particle P to an observer in the moving 
frame. 


(d) 



Angular acceleration of the moving frame to an observer in 
the inertial frame. 


(e) 


2 5 x 


Vi 

dt 


Coriolis acceleration of P generated due to rotation. 


(f) nx (nxr) Centripetal acceleration of P when fixed in position in the 

moving frame, and the moving frame is rotating at constant 
angular velocity. 


For an observer in the moving frame (with no translation) the terms (a), (b) 
and (d) in equation (A3. 9) would disappear and equation (A3. 9), in moving 
frame becomes: 


d l - d l^l - - - 

x + J 2 x(s 2 xr.) s 0 

dr dt 1 


(A3. 10) 
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writing 


dj 2 rj 

j- * A, acceleration of particle P 

dt 


and 


d i p i 

dt 


7, velocity of particle P 


in the rotating frame equation (A3. 10) becomes: 

A + 25 x 7 + n x (n x fj) =0 


(A3. 11) 


So an observer in the rotating frame can use the components of the 

he ,« es <;> f° r °v/Dt, Dw/Dt On momentum equation), 

provided he includes the extra apparent body forces of per unit mass, 

-2 [n x 7 + n x (n x rj) ] 


where 


= + Jfiy + V = 1u + jv + kw, rj = ix + jy + kz 

9 T r,te ? due f° ****«■ rotation in momentum equations 
10 [2^ x t 7 e + t ^ rmS (^ h ° Wn *jj 0Ve shou ^ multiplied by density (-), 

These terms go into generalized equation in TEACH as source terms (Sj^j) , i.e. 


^ » -p[2S x 7 + n x (jj x ? x ) 


(A3. 12) 


Resolving equation (A3. 12) in 1, j, and k directions gives: 

p2(n x v) = -2ip[ws^ - vn 2 ] + 2jp[wj^ - ufi 2 3 - 2kp[vfi x -ufi y ] (A3. 13) 

-p£ 2 x (n x r 2 ) = -ipfoy (yn x - xj2 y) - fi z (xn z -zn x ) ] 

+ jp [n x (ys 2 X - XI 2y) - 12 z (zn y - yn z ) ] 

- kp E « x (XS2 Z - ZI2 X ) -I2y (Zl2 y - yi 2 z ) ] (A3. 14) 

Adding (A3. 13) and (A3. 14) we get: 

= ~~ p ^ x 7 + n x (S x rj) ] 
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= i [ p(v - ZS^ + xn z ) J 2 Z - p(w - Xj2y + yn x )n y + pVfi z pwn y ] 

+ j [ p(w - xj^ + yn x ) n x - p(u - yn z + zs^)q z + pws^ - pus^ ] 

+ k [ p(u - m z + zfiy) j2 y - p(v - zj^ + xn z )a x + pun^ - pvj^ ] 

(A3. 15) 

Contribution of Rotation in Stagnation Enthalpy Equation (H rot) 

Contribution of rotation to the stagnation in enthalpy can be found by forming 
dot product of the velocity vector ( v) and the extra apparent body forces 
generated by the system rotation. This becomes: 

H rot = pV ' [ - (2 fiV + ft x (ft x r^))] 

= P V • [ - (fi x (fix rj))] 

= - P V ' [ft (ft ’ ? x ) - r 1 (ft • ft) ] 

These extra terms generated due to rotation go into the generalized TEACH 
equations as extra source terms (S h) » l.e. 

S^ H = - pV [ fi (ft ’ r x ) - r x (ft * ft) ] 

s nH = - ’ c «(fi x x + + fi z 2 > - + + « z 2 ) ] 

= - p (iu + jv + kw) ' [ (ifi x + jfiy + kn z ) (fi x x + + fi 2 z) 

- ( lx + jy + kz) (fi 2 + fi 2 + fi 2 ) ] 

a y 4. 

= - p C ( ufi x + vfi y + wfi^ ( f^x + fi^y + fi 2 z) 

- (ux + vy + wz)(fi 2 + fi 2 + fi 2 ) ] 

a y z 

collecting the u, v, w terms and averaging gives: 

2 2 2 2 
SfiH = - pu [ fi x x + n x ^y + fi A z - n x x -n y x -fi 2 x ] 
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- pv [n^x + n y 2 y + nn z z - m z - y^ 2 - y^ 2 ] 

- pw [ n^x + + n z 2 z - z^ 2 - zj^ 2 - zn z 2 ] 

or 

s h = - C fiy - n y x) + n z (n x z - n z x) ] 

- pv [ S2 x (n^x - + n z (n^z - n^) ] 

- pw [ n x (n z x - j2 x z) + n y (n^ - n y z) ] 

(A3. 16) 


The expression In (A3. 15) has been modeled Into 3D-TEACH separately In u, v 
and w momentum equations, l.e., 1-component In u-equatlon, j-component In 
u-equation, and k-component In w-equation. 

The expression in (A3. 16) has been modeled in 3D-TEACH code as shown above. 
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